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I .   INTRODUCTION 

If  F(t)  is  an  integrable  function  of  the  real  variable 
t,  the  function  f(s)  defined  by  the  integral 


/-st 
e 


f(s)  =  /    e    F(t)dt  (1) 


is  called  the  Laplace  transformation,  or  the  Laplace  trans- 
form, of  F(t)  and  is  frequently  indicated  by  the  notation 


f(s)  =£[F(t)]  (2) 


The  Laplace  transformation  is  a  linear  integral  trans- 
formation which  is  widely  used  in  applied  mathematics  and 
technology.   A  typical  application  is  as  follows.   An  un- 
known function  F(t)  satisfies  a  certain  differential  equa- 
tion and  specified  end  conditions.   Employing  theorems 
applicable  to  the  Laplace  transform,  a  function  f(s)  is 
determined  as  the  solution  of  an  algebraic  equation  which 
corresponds,  in  an  appropriate  fashion,  to  the  differential 
equation  and  end  conditions  of  F(t).   Then  the  problem  is 
reduced  to  the  following  task:   having  established  f(s), 
it  is  necessary  to  determine  the  function  F(t)  for  which 
f(s)  =£[F(t)]. 


This  problem  is  called  finding  the  inverse  (Laplace) 
transform  of  the  function  f(s).   This  is  customarily 
written  as 

F(t)  =<£_1[f(s)]  (3) 

It  may  be  shown  that  the  following  inversion  integral 


F(t)  =  2?I  /  eSt  f(s)ds  (4) 


accomplishes  the  desired  inversion.   In  this  expression 
s  is  the  complex  variable 

s  =  x  +  iy  (5) 

and  the  integration  is  along  a  suitable  path  in  the  complex 
s-plane. 

If  in  fact  there  is  a  function  F(t)  such  that  the  given 
f(s)  =cL[F(t)],  then  the  function  determined  by  equation 

(4)  is,  indeed,  the  function  F(t).   The  result  is  unique 
in  the  following  sense.   If 


e"St  F1(t)dt  (6) 


and 


F2(t)  "  2Ti  /  ^  f(s)ds  (7) 


then  F,(t)  and  F„(t)  differ  only  on  a  set  of  Lebesgue 
measure  zero. 

In  this  thesis  numerical  methods  will  be  utilized  to 
perform  the  integration  of  equation  (4)  and  functions  for 
which  the  inversion  is  successful  will  correspond  to 
functions  F(t)  which  are  at  least  piecewise  continuous, 
having  only  isolated  discontinuities,  if  any  are  present 
at  all.   If  F(t)  is  discontinuous  for  t  =  a,  then  the 
recovered  function  satisfies  the  condition 

F(a)  =  -|  Lim[F(a+e)  +  F(a-e)]  (8) 

e  -»■  0 

The  path  of  integration,  indicated  by  C  in  equation 
(4)  is  frequently  described  as  an  infinitely  long,  vertical 
straight  line  x  =  constant,  to  the  right  of  any  singulari- 
ties of  f(s).   (A  theorem  of  complex  variables  shows  that 
f(s)  is  analytic  in  its  region  of  convergence,  a  right-hand 
plane,  except  for  isolated  singularities.)   This  straight 
vertical  contour  is  frequently  called  the  Bromwich  contour. 

Usually  inversion  of  a  Laplace  transform,  that  is  find- 
ing F(t)  =  oC~    [f(s)],  is  accomplished  by  use  of  a  table  of 
transform  pairs.   One  of  the  most  complete  of  these  is  that 
of  Roberts  and  Kaufman  [14] .   Other  such  tables  are  listed 
as  references  [5],  [6],  [7],  and  [13], 

Several  theorems  concerning  the  Laplace  transform  extend 
the  usefulness  of  such  transformation  pair  tables.   For 
example 
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<£  1[f1(s)  +  f2(s)]  =£   1[f1(s)]  +£   1[f2(s)] 

(9) 

However,  it  is  not  an  infrequent  occurrence  that  one  is 
called  upon  to  find  F(t)  in  a  case  where  the  tables  and 
the  theorems  are  of  no  direct  help.   There  has  developed 
a  considerable  literature  concerning  this  problem. 

One  obvious  method  is  to  attempt  to  expand  the  given 
f(s)  as  a  series  of  functions  t(s)for  which  the  inverses 
are  known.   Then  using  the  linear  property  indicated  in 
equation  (9),  the  result  may  be  expressed  as 

c£"1[f(s)]  =dC"1[Zaifi(s)]  =  Zaj«t'1[fi(s)]         (10) 

The  integration  indicated  in  equation  (4)  can  also  be 
performed  by  various  analytical  and  approximate  methods. 
In  particular,  one  can  perform  the  integration  along  the 
Bromwich  contour  by  a  numerical  procedure.   Much  of  the 
literature  is  devoted  to  such  methods. 

In  general  these  methods  amount  to  the  following.   A 
finite  number  of  points  along  the  Bromwich  contour  are 
determined  according  to  some  law  and  the  integrand  is 
evaluated  at  such  points.   There  are  associated  weight  fac- 
tors and  the  integral  is  approximated  as  a  sum  of  products 
of  integrand  values  times  weight  factors.   One  of  the  most 
widely  used  of  such  methods  is  described  by  Salzer  [15] , 
[16],  and  [17].   Other  procedures  of  the  type  are  treated 
in  [18],  [20],  [2],  [8],  [9],  [10],  and  [21]. 
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It  is  interesting  to  note  that  the  same  classical  poly- 
nomials which  are  useful  in  the  expansions  indicated  by 
equation  (10)  are  also  encountered  in  locating  points  for 
evaluation  using  a  numerical  integration  scheme  so  that 
the  final  algorithm  may  be  the  same  even  though  the  original 
motivation  was  quite  different. 

Other  methods  which  do  not  fall  clearly  into  the 
classes  of  methods  discussed  above  have  been  presented  by 
Widder  [22],  and  Bellman,  Kalaba,  and  Lockett  [3]. 

These  approximate  methods  have  not  always  been  success- 
ful.  Hiep  [12]  employed  Salzer's  method  for  numerical  in- 
version on  a  problem  of  heat  transfer  in  porous  media  and 
found  it  difficult  to  obtain  results  which  did  not  exhibit 
physically  incorrect  behavior,  e.g.,  in  the  vicinity  of  a 
point  of  engineering  interest,  the  temperature  of  the 
medium  fell  below  sink  temperature  and  rose  above  source 
temperature.   In  similar  work  on  a  conjugated  heat  transfer 
problem  Zargary  [23]  encountered  the  same  difficulty,  and, 
after  devoting  considerable  time  and  attention  to  the  prob- 
lem, was  forced  to  abandon  the  use  of  numerical  inversion. 

To  overcome  such  difficulties  in  problems  of  these 
kinds,  we  have  reverted  to  direct  numerical  integration  of 
equation  (4)  along  appropriate  contours  in  the  s-plane. 
Our  procedures,  which  appear  thus  far  to  be  efficient  and 
reliable,  are  described  in  what  follows. 
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II.   DISTORTED  CONTOURS 

We  have  found  that  great  advantage  accrues  from  distort- 
ing the  Bromwich  contour  in  ways  which  will  subsequently 
be  described. 

According  to  theorems  concerning  analytic  functions  of 
a  complex  variable,  the  Bromwich  contour  can  be  distorted 
arbitrarily  as  long  as  the  process  of  distortion  does  not 
cause  the  contour  C  to  cross  over  a  singularity  of  f(s)  and 
as  long  as  the  two  ends  of  C  extend  to  infinity  in  a  manner 
which  preserves  the  convergence  of  the  integral  in  equation 
(4). 

Thus,  with  the  supposition  that  f(s)  is  analytic,  with 
the  exception  of  isolated  singularities,  and  is  real  valued 
for  real  s,  let  g(s)  represent  the  integrand 


g(s)  =  eSt  f(s)  =  u(s)  +  iv(s)  (11) 


where  u(s)  and  v(s)  are  the  real  and  imaginary  harmonic 
components  of  g(s).   Furthermore,  assume  that  the  contour  C 
consists  of  the  two  symmetrical  parts  C   and  C    as  shown 
in  Figure  1.   Consideration  of  the  contributions  tooL   [f(s)j 
over  arcs  ds-.  on  C   and  ds2  on  C   yields  the  following: 


dx2 

= 

-dx. 

dy2 

= 

dyx 

U2 

~ 

ul 

v2 

— 

■vl 

(12) 
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-> REAL  AXIS 


FIGURE  1.   Distorted  Bromwich  Contour.   Contour  is 
distorted  into  C*  and  its  reflection  C  * 


Thus, 


27i^ldsl  +  g2ds2] 


¥[v1dx1  +  Uldyi] 


(13) 


and,  therefore, 


F(t)  =  - 

TT 


[vdx  +  udy] 


(14) 


Equation  (14)  is,  of  course,  further  simplified  if  C 
is  selected  to  be  a  path  along  which  v=0.   Additional 
theoretical  reasons  exist  for  such  a  selection.   Such  a  path 
will  subsequently  be  referred  to  as  a  path  of  steepest 
descent.    Numerical  integration  along  such  a  path,  or  a 
close  facsimile  thereof,  was  done  by  Esch  [11]  in  1957  and 
Carrier,  Krook,  and  Pearson  [4]  recommended  further 
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exploration  of  this  procedure.   However,  no  further  such  in- 
vestigation seems  to  have  been  done  until  the  present  thesis, 
In  Chapter  V,  a  description  of  an  algorithm  for  performing 
numerical  integration  along  a  path  C  ,  for  which  v=0  over 
a  substantial  portion  of  its  range,  will  be  presented. 
Almost  by  accidental  discovery,  the  present  work  has 

revealed  that  it  is  computationally  more  efficient  simply 

* 
to  make  an  arbitrary  choice  for  the  contour  C  which  satis- 
fies the  following  criteria: 

(1)  C   is  a  smooth  continuous  curve  of  the  form 

s(p)  =  x(p)  +  iy(p)  (15) 

where  p  is  a  real  parameter  0   p   °° . 

(2)  C   lies  to  the  right  of  all  singularities  of  f(s). 

(3)  C   does  not  extend  sufficiently  far  to  the  right 

to  cause  computational  difficulties  due  to  the 

,       st 
factor  e 

(4)  C   approaches  infinity  in  such  a  way  that  x 
approaches  -°°  and  y  approaches  °°. 

(5)  Oscillation  of  the  integrand  is  not  excessive. 
(Otherwise  there  may  be  loss  of  accuracy  due  to 
the  sampling  rate  of  the  integration  algorithm 
seriously  mismatching  variations  in  the  integrand 
and/or  positive  and  negative  contributions  nearly 
cancelling  each  other.) 

Such  a  contour  will  be  referred  to  as  a  simple  parameterized 
contour.   Discussion  of  how  to  assure  satisfying  criteria 
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(2),  (3),  and  (5)  will  be  deferred  until  after  the  develop- 
ment of  the  actual  algorithm. 

The  integral  shown  in  equation  (14)  is  approximated  by 
the  simple  trapezoidal  sum 

M 

F<t}  =  2T  Z  {[v(sk)+v(sk+i)]txk+i  "  xlJ 
k=l 

+  [u(sk)+u(sk+1)][yk+1-  yk]}         (16) 

where 

Sk  =  Xk  +  iyk  =  x(pk)+  iy(pk}        (17) 

The  arbitrarily  chosen  functions  x(p)  and  y(p)  are  monotoni- 
cally  increasing  and  p,  ,  p2,««.,  P   is  an  increasing 
sequence . 

A  very  simple  example,  but  one  which  has  been  employed 
very  successfully,  is  given  by  the  following  equations 

x  =  A  -  Bp2 

y  =  P  (18) 

p  =  0,  A,  2A, .. . 

The  real  number  A  is  selected  so  as  to  satisfy  criteria  (2) 
and  (3)  above.   The  real  number  B  has,  usually,  been 
assigned  a  value  of  one. 

The  sum  includes  only  a  finite  number  of  terms.   It  is 
necessary  to  provide  an  appropriate  criterion  for  terminat- 
ing the  process  of  summation.   This  will  be  discussed  later 
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In  order  to  monitor  the  possibility  of  inaccurate 
evaluations  resulting  from  oscillations,  in  addition  to  the 
sum  F(t)  given  by  equation  (16),  the  sum  G(t)  given  by 


G(t)  = 


_1_ 

2tt 


k=l 


([v(sk+1)-v(sk)][xk+1-xk] 


+  [u(sk+1)+u(sk)][yk+1-yk]} 


(19) 


is  also  obtained.   This  is  the  sum  of  the  absolute  values 
of  the  addends  to  the  sum  (16).   Additionally,  the  number 
of  times  successive  addends  are  of  opposite  signs  is 
recorded. 

The  present  numerical  experiments,  to  be  described 
later,  were  generally  so  successful  that  the  effect  of 
oscillation  could  not  be  seen.   However,  in  some  cases 
where  poor  contours  C   were  deliberately  chosen,  it  was 
found  that  inaccurate  results  were  obtained  if  G(t)  was 
many  orders  of  magnitude  (e.g.,  10  )  times  as  great  as  F(t) 
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III.   THE  COMPUTER  PROGRAM  AND 
ITS  PRINCIPAL  SUBROUTINES 

Previous  discussion  proposed  two  alternatives  in  terms 
of  advantageous  distortion  of  the  Bromwich  contour  upon 
which  to  perform  the  numerical  integration  of  equation  (16). 
The  simple  parameterized  curve  has  been  experimentally  found 
to  be  the  more  efficient  of  the  two  methods  and  its  algo- 
rithm will  be  developed  at  this  time.   The  steepest  descent 
procedure  will  be  treated  in  Chapter  V. 

The  computer  program  and  subprograms  which  implement 
the  simple  procedure  developed  in  the  preceding  chapter  are 
all  written  in  the  FORTRAN  language,  using  double  precision 
arithmetic.   They  have  been  tested  and  debugged  on  an  IBM 
360/67.   For  machines  having  a  larger  mantissa,  single  pre- 
cision arithmetic  might  prove  to  be  satisfactory. 

SUBROUTINE  VALUE  is  a  user  supplied  routine  which  imple- 
ments computation  of  the  g(s)  given  by  equation  (11)  for  the 
desired  transform  f(s).  This  subroutine  calculates  the  u(s) 
and  v(s)  corresponding  to  a  given  x  and  y. 

Either  complex  or  real  arithmetic  may  be  utilized  in  this 
program,  although  complex  arithmetic  offers  a  decided  advan- 
tage in  convenience  and  simplicity.   However,  the  user 
should  be  aware  of  the  potential  dangers  in  unintentional 
exchange  of  sheets  of  a  Riemann  surface  when  investigating 
a  multi-valued  function   More  will  be  said  about  this  later. 
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Also,  although  FORTRAN  manuals  state  or  imply  that  the 
FUNCTIONS  DREAL  and  DIMAG  are  implemented  in  all  compilers, 
this  may  not  actually  be  the  case.   Accordingly  it  is  prudent 
to  add  the  following  statements  after  SUBROUTINE  VALUE. 


FUNCTION  DREAL (CPLX) 

REAL* 8  CPLX(2) 

DREAL  =  CPLX(l) 

RETURN 

END 

FUNCTION  DIMAG (CPLX) 

REAL* 8  CPLX(2) 

DIMAG  =  CPLX(2) 

RETURN 

END 


SUBROUTINE  CURVE  is  also  a  user  supplied  routine  which 
calculates  the  x  and  y  values  for  each  consecutive  incrementa- 
tion of  the  parameter  p  of  equation  (17),  thus  locating  a 
point  on  the  simple  parameterized  contour.   Equation  (18) 
provides  an  illustrative  and  highly  effective  example  of 
parabolic  form  which  has  been  employed  almost  exclusively 
during  this  investigative  work. 

SUBROUTINE  INCREP  formulates  the  increasing  sequence  of 
the  parameter  p  as  given  by  the  following  equation 

p  =  p  +  A  (20) 

Limited  analysis  with  other  than  equally  spaced  incrementa- 
tion of  points  along  the  contour  of  integration  did  not  ex- 
hibit any  enhancing  capabilities  and  was  not  continued. 

The  main  program  accomplishes  the  integration  of  equation 
(.14)  using  a  trapezoidal  summation  as  described  by  equation 
(16).   In  addition,  the  absolute  value  of  the  integrand  is 
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summed  in  accordance  with  equation  (19).   Also,  the  number 
of  times  a  sign  changes  occurs  between  successive  evaluations 
of  the  addend  of  equation  (16)  is  recorded.   These  features 
afford  the  user  an  opportunity  to  monitor  potential  oscilla- 
tion effects  upon  the  accuracy  of  the  results.   Clock  time 
and  the  number  of  calls  to  SUBROUTINE  VALUE  required  for 
performing  the  numerical  contour  integration  are  maintained 
as  measures  of  the  computational  efficiency  of  the  numerical 
inversion  procedure. 

Termination  of  numerical  integration  occurs  when  each 
of  N  successive  evaluations  of  the  addend  of  equation  (19) 
has  magnitude  less  than  a  specified  epsilon.   We  have  usually 
used  N  =  5.   This  is  a  very  stringent  requirement  and  could 
probably  be  relaxed  resulting  in  some  saving  of  computational 
time  and  with  but  minor  loss  of  accuracy. 

In  addition  to  the  termination  criterion  described  in 
the  previous  paragraph,  other  program  inputs  include  the 
increment  A  of  equation  (18),  appropriate  numerical  values 
of  any  constants  of  the  function  f(s),  the  starting  position 
A  of  equation  (18)  and  the  value  of  t. 

The  program  output  prints  the  starting  position  of  the 
integration  contour,  the  values  given  by  equations  (16)  and 
(19),  the  final  values  of  x  and  y  at  the  termination  of  inte- 
gration, the  number  of  changes  of  sign  between  successive 
addends  in  equation  (16),  the  total  number  of  calls  made  to 
SUBROUTINE  VALUE  during  the  process,  and  the  clock  time  re- 
quired for  numerical  integration. 


20 


As  an  illustration  consider  the  following  example  where 


f(s)  =   2  -   2  (21) 

s   +  ct 


A  suitable  form  for  SUBROUTINE  VALUE  is  as  follows: 

SUBROUTINE  VALUE  (X,  Y,  T,  AL,  U,  V,  MMMMM) 
C      THIS  IS  TABLE  ENTRY  NUMBER  09  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 
IMPLICIT  REAL*8  ( A,B,D-H ,0-Z) 
IMPLICIT  COMPLEX*16  (C) 
MMMMM  =  MMMMM  +  1 
ZERO=0.D+0 

CS=DCMPLX(X,Y)  ,„  .     _. 

CAL=DCMPLX(AL,ZERO)        ^ote  [      The  parameter  a  is 
CT=DCMPLX(T,ZERO)  denoted  by  the  FORTRAN 

CDEN=CS**2+CAL**2  variable  AL. ) 

C=CS/CDEN 
CEXP=CDEXP ( CS  *CT ) 
C=C*CEXP 
U=DREAL(C) 
V=DIMAG(C) 
RETURN 
END 

If  the  contour  indicated  by  equation  (18)  were  chosen  and 

numerical  inversion  were  performed  upon  the  f(s)  given  by 

equation  (21)  for  the  case  where  A  =  0.3,  A  =  0.1,  a  =  0.125, 

N  =  5  and  e    =   l.D-11,  the  evaluation  of  F(2tt)  would  yield  the 

following  output 

AA    F(2tt)        |G(2tt)|      Xf      Yp    LOSC    MMMMM 
0.3   7.07107D-1   1.32635D+0   -5.46   2.4    4      25 

where  AA,  LOSC,  and  MMMMM  are  the  FORTRAN  program  names  for 
the  contour  starting  position,  the  number  of  sign  changes 
between  successive  addends,  and  the  total  number  of  calls 
made  to  SUBROUTINE  VALUE,  respectively.   Similar  results 
for  the  almost  100  functions  f(s)  tested  during  the  course 
of  this  research  are  listed  in  Appendix  B. 
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IV.   AUXILIARY  SUBROUTINES 

Chapter  III  briefly  mentioned  the  potential  dangers  in- 
herent in  the  utilization  of  complex  arithmetic  within  SUB- 
ROUTINE VALUE  when  the  function  f(s)  is  multi-valued. 

If  f(s)  has  branch  points  or  winding  points  it  is  impor- 
tant to  assure  that  one  remains  on  the  same  branch  of  the 
function  as  one  proceeds  along  the  path  of  integration. 
Available  FORTRAN  operations  do  not  assure  this  and  it  is 
necessary  to  employ  specially  programmed  algorithms.   The 
present  discussion  is  limited  to  the  two  simple  cases,  rais- 
ing of  a  complex  number  to  a  non-integer  power  and  taking 
the  logarithm  of  a  complex  number. 

The  essence  of  the  matter  lies  in  the  FORTRAN  operation 
DATAN2(y,x),  or  its  equivalents,  which  must  be  relied  upon 
to  provide  the  argument  of  a  complex  number  in  polar  form. 
This  operation  always  provides  a  result  in  the  range 
-7T  <  9  <_  7T ,  so  that  if  Arg  (z)  actually  passes  through  the 
values  (2n+l  )tt,  for  integer  n,  the  FORTRAN  produced  value 
of  Arg(z)  experiences  a  jump  of  ±  2tt  .   If  one  is  determining 
the  logarithm  of  z  or  a  non-integer  power  of  z,  shifting  of 
branch  will  take  place  unless  the  continuity  of  Arg(z)  is 
restored. 

SUBROUTINES  CPOWER  and  CLOG  provide  for  this  continuity 
by  monitoring  changes  in  Arg(z)  and  adding  or  subtracting 
2niT ,  whenever  it  is  appropriate  to  do  so,  to  the  value  of 
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DATAN2(y,x).      As   they    are   presenting  written,    each  provides 
for  up   to    five   separate   calls,    within   SUBROUTINE   VALUE,    to 
either   SUBROUTINE    CPOWER   or   SUBROUTINE   CLOG. 
These   subroutines   are    listed  below, 

SUBROUTINE    COQWER    ( 11 N, ZCUT, P, J , N CALL ) 

IMPLi:iT    REAL*8    (A-H,0-Z) 

CQMPLEX*16     ZIN,     ZOUT 

DIMENSI CN    T0LD(5) ,    NCALL(5),    MPK5) 

IF(MCALL(J) .E0.1)    NPI(J)=0 

OH   =    0  .  D+0 

X    =    DREAL( ZIN) 

Y  =    DIM AGfZIM ) 

R    =    DSORT(  (X*X)     +    (Y*Y  i  ) 

TRQVJ3     =    DATAN2(Y,X) 

IF(NCALL( J  )  .E0.1)    GO    TO    5 

PROD    =    ?RONG*TOLD(JJ 

IF(  PPOD.GT.OH.OR.  X.GT.CH)     GO    TO    5 

IF(TRONG.LT .OH)    MPI(J)    =    NPI(J)    +    2 

IF{TRONG.G".OH)     NPI(J)     =    NPI(J)    -    2 

5     RP     =      R**D 

EN   =    MPI  (J  ) 

TRITE    =    TRONG    +     ( EN*3.14159265358979324D+0) 

TP    =     TRITE*P 

X    =    RP*CCOS  (TP) 

Y  =    RP*DSIM(TP) 
ZOUT    =    DCMPLX(X,Y) 
T0LD(J)     =    TCPNG 
NCALL(  J)     =    MCALLt J)     +    1 
RETURN 

END 


SUBROUTINE    CLOG    (ZIN,     Z3UT,     J,    ^CALL) 

IMPLICIT    REALMS     (A-H,Q-Z) 

C3M?LEX*16    ZIN,    ZOUT 

DIMENSION    T0LD(5),    MCALL(5),     NPK5) 

IF(MCALL(J)  .EQ.1J    NPKJ)    =0 

^H   =     O.D+O 

X    =    DREAL(ZIN) 

Y   =    DIMAG( ZIN) 

R    =    DS0RT( { X*X)     +    { Y*Y) ) 

TRONG    =    DA-*- AN2  (  Y  ,  X  ) 

IF(  M:^LL(  J)  .EQ.  1)     GO    TO    5 

PROD    =    TRONG*TOLD(J ) 

IF(PROD.GT.OH.OR.X.GT.OH)    GO    to    5 

IF(TRONG.LT.OH)     NPKJ)    =   NPI(J)    +    2 

IF(TRONG.GT.OH)    NPKJ)    =    NPKJ)    -    2 

REAL    =    0L0G(R) 

EN    =    M P  I  (  J  ) 

TRKE    =    TRONG    +    ( EN*3.14159265353979324D+0) 

ZCUT    =     DCMPLX(R£AL, TRI^E) 

TOLD( J)    =    trong 

MCALL (J )    =    MCALL( J)     +    1 

RETURN 

END 
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If  either  algorithm  is  to  be  utilized  by  SUBROUTINE  VALUE, 
the  first  requirement  is  to  zero  the  applicable  integer 
array,  NCALL  or  MCALL,  within  the  main  program.   Then,  in 
SUBROUTINE  VALUE,  the  calls  to  SUBROUTINE  CPOWER  are  of 
the  form 

CALL  CPOWER ( Z I N , ZOUT , P ,  J,  NCALL ) 

where  ZIN  is  the  C0MPLEX*16  operand,  ZOUT  is  the  C0MPLEX*16 
output,  P  is  the  REAL*8  power,  and  J  is  the  INTEGER*4  index 
which  tells  whether  this  is  the  first,  second,...,  call  to 
SUBROUTINE  CPOWER  within  SUBROUTINE  VALUE.   The  relationship 
obtained  is 

ZOUT  =  (ZIN)P  (22) 

SUBROUTINE  CLOG  is  used  similarly.   Calls  made  from  SUBROUTINE 
VALUE  are  of  the  form 

CALL  CLOG (ZIN, ZOUT, J, MCALL) 

where 

ZOUT  =  Ln(ZIN)  (23) 

J  plays  the  same  role  as  in  SUBROUTINE  CPOWER.   There  is  no 
P.   It  should  be  noted  that  the  array  MCALL  replaces  NCALL. 
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V.   THE  STEEPEST  DESCENT  CONTOUR 

The  theoretical  advantages  of  a  steepest  descent  contour 
have  been  mentioned  oreviously.   Consider  a  contour  in  the 
complex  plane  along  which  v=constant.    As  shown  in  Figure 
2,  let  the  |  axis  be  tangent  to  the  contour  with  the  n  axis 
perpendicular  to  the  £  axis  in  the  same  sense  that  y  is  perpen- 
dicular to  x. 


■>  REALAX.I5 


Figure  2.   Steepest  Descent  Contour.   Contour  along 

which  v  is  constant.   Axis  £  is  tangent  at 
point  P.   Axis  n  is  normal. 


Recall  equation  (11), 


g(s)  =  est  f(s)  =  u(s)  +  iv(s)  (11) 


By  the  Cauchy-Riemann  relationships 


8u   _3v      9v    3u 
Along  a  v  =  constant  contour 


(24) 


|V  =  o  ;   |H  =  0  (25) 
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Such  a  v  =  c  =  constant  contour  divides  the  plane,  locally 
at  least,  into  two  regions:   R,  with  v  >  c  and  R~  with  v  <  c. 
Consequently  along  such  a  contour  -r—  is  one  signed,  and, 
therefore,  «  is  one  signed.   u  •*■  0  along  the  v  =  c  contour, 
in  the  direction  of  positive  E, .      Therefore,  there  is  no 
oscillation  in  u  and,  obviously  also,  there  is  none  in  v. 
Since  u  -*•  0  and  v  •*■  0  along  a  suitable  integration  contour, 
the  selection  of  c  =  0  assures  that  this  results. 

This  second  program  to  be  discussed  does  partially  follow 
a  path  of  steepest  descent.   First  it  integrates  vertically 
upward  from  a  chosen  point  on  the  real  axis  (along  a  Bromwich 
contour)  until  intersection  with  a  specified  v  =  0  contour 
occurs.   Then  it  employs  a  tracking  algorithm  to  obtain  suc- 
cessive points  on  this  contour  along  which  it  then  integrates, 
moving  to  the  left. 

SUBROUTINE  VALUE  has  the  same  role  as  that  described  in 
Chapter  III,  namely,  it  calculates  the  u(s)  and  v(s)  in 
accordance  with  equation  (11)  corresponding  to  a  given  x  and 
y.   For  the  initial  part  of  the  algorithm,  the  integration 
along  the  Bromwich  contour,  the  value  of  x  is  constant,  and 
y  is  increased  in  increments  by  Ay.   Numerical  integration 
in  accordance  with  equation  (16)  is  conducted  just  as  it  was 
done  in  the  previous  case. 

On  the  real  axis,  obviously,  v  =  0.   With  an  initial 
incrementation  in  the  vertical  direction,  v  assumes  a  value 
other  than  zero.   The  sign  of  v  is  recorded  and  the  incrementa- 
tion process  is  continued  until  the  sign  of  v  changes.   This 
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constitutes  a  contour  crossing.   The  procedure  continues 
until  a  predetermined  number  of  such  crossings  are  obtained, 
and  then  a  simple  linear  interpolation  is  performed  to  locate 
the  specific  ordinate  value  y  for  which  v  =  0.   Usually  one 
might  wish  to  use  the  first  v  =  0  contour  above  the  real 
axis,  but  the  second,  third,  or  other  such  contour  may  be 
selected. 

The  last  (y.,v.)  values  prior  to  the  desired  contour 
crossing  and  the  first  (y.,v.)  values  after  crossing  are 
utilized  in  the  linear  interpolation  process  to  predict  YNEW, 
for  which  v  =  0.   SUBROUTINE  VALUE  provides  the  actual  VNEW 
value  corresponding  to  YNEW.   Depending  upon  the  sign  of 
VNEW,  the  appropriate  coordinate  pair  (y*>v±)    or  (y->v-)  is 
replaced  by  (YNEW, VNEW)  and  the  interpolation  is  repeated 
with  this  refinement.   Termination  of  the  interpolation  occurs 
when  YNEW  results  in  a  magnitude  of  VNEW  which  is  less  in 
magnitude  than  a  prescribed  tolerance  limit. 

Once  the  values  of  x,y,u,  and  v  have  been  established 
for  the  first  point  on  the  desired  v  =  0  contour  in  the  pre- 
vious process,  the  second,  or  tracking  portion  of  the 
algorithm  commences. 

SUBROUTINE  XMARCH  is  called  to  provide  a  new  coordinate 
pair  (x,y)  in  the  negative  sense  of  the  real  axis  from  the 
known,  or  old,  coordinate  pair  (x  ,y  ).   This  is  accomplished 
by  the  following  relationships 

x  =  xQ  -  R  cos(6)  (26) 

y  =  yo  +  R  sin(9)  (27) 
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where  (R,9)  are  the  polar  coordinates  of  a  search  pattern 
measured  from  the  known  coordinates  (x  , y  )  on  the  v  =  0 
contour.   The  angle  9  is  set  at  zero  radians  for  the  first 
call  to  SUBROUTINE  XMARCH.   SUBROUTINE  VALUE  is  then  called 
to  provide  the  value  of  v  at  this  point.   The  sign  of  v 
establishes  upon  which  side  of  the  v  =  0  contour  the  new  point 
lies. 

SUBROUTINE  ROTATE  is  then  called  to  increment  the  angle 
theta  in  an  appropriate  sense.   A  sweeping  search  is  then 
conducted  by  SUBROUTINES  ROTATE,  XMARCH,  AND  VALUE  until 
a  sign  change  is  detected  for  v  values  between  increments, 
thus,  indicating  a  crossing  of  the  v  =  0  contour. 

The  last  (9j,v.)  coordinate  pair  prior  to  a  crossing  and 
the  first  (9., v.)  coordinate  pair  after  a  crossing  are 
utilized  in  a  linear  interpolation  scheme,  analogous  to  that 
of  the  first  program  segment,  to  refine  the  location  of  the 
new  point  on  the  v  =  0  contour. 

Once  the  second  point  on  the  v  =  0  contour  is  established, 
the  tracking  algorithm  may  be  enhanced  by  predicting  the 
third  point  on  the  contour  to  lie  at  the  same  angle  9  from 
the  second  point.   As  before,  the  sign  of  v  at  this  new  loca- 
tion permits  a  rotational  search  to  be  conducted  in  the 
appropriate  sense.   For  small  values  of  the  radius  vector 

R  (R=.025  is  built  into  the  program,  but  may  be  changed  by  the 
user),  this  approximation  has  been  found,  in  general,  to  be  a 
reasonable  one. 

The  addends  of  equations  (16)  and  19)  are  evaluated  at 
each  step  in  the  marching  process  within  the  main  program 
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and,  once  again,  termination  of  the  numerical  integration 
occurs  when  each  of  N  successive  addend  contributions  is 
less,  than  a  prescribed  epsilon. 

Required  inputs  to  the  main  program  are  the  index  speci- 
fying the  desired  v  =  0  contour,  the  starting  position  on 
the  real  axis,  the  increment  of  the  vertical  march  along  the 
Bromwich  contour,  the  termination  criterion,  any  constants 
associated  with  the  function  f(s),  and  the  value  of  t.   These 
inputs  are  represented  by  the  program  variables  NREG,  AA, 
DELY,  EPS,  (AL,BE,...),  and  T,  respectively. 

The  simple  parameterized  contour  offers  one  decided  ad- 
vantage over  the  steepest  descent  contour:  each  call  made  to 
SUBROUTINE  VALUE  is  productive,  in  the  former  case,  in  the 
sense  that  it  is  utilized  to  compute  a  contribution  of  the 
addend  of  equation  (16).   The  latter  algorithm  requires  calls 
to  SUBROUTINE  VALUE  to  perform  the  linear  interpolation  pro- 
cedures, which  are  non-productive  in  terms  of  addend  computa- 
tion.  This  represents  increased  computer  time.   The  present 
investigation  has  shown  that  in  each  of  the  numerous  cases 
which  have  been  examined,  the  simple  parabolic  path  indicated 
by  equation  (18)  is  more  efficient  and  fully  as  accurate  as 
a  steepest  descent  contour. 
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VI.   SHANKS '  ACCELERATOR 

The  effectiveness  of  a  family  of  non-linear  sequence  to 
sequence  transformations  in  accelerating  the  convergence  of 
(some)  slowly  convergent  sequences  and  in  inducing  conver- 
gence in  (some)  divergent  sequences  has  been  reported  by 
Shanks  [19].   If  (A  }  (n=0, 1 , 2 , . . . )  is  a  sequence  of  numbers, 

we  form  a  sequence  of  sequences  (A    }  which,  for  convenience, 

k  ,  n 

we  write  as  (A(k,n)}.   The  integer  k  indicates  the  "order" 
of  the  sequence,  with  A   =  A(0,n);  the  integer  n  indicates 
the  position  of  the  term  in  the  sequence.   The  rule  for  form- 
ing the  sequence  (A(k+l,n)}  from  the  sequence  (A(k,n)}  is 

A(i,j)  =0   if   j  <  2i 

A(k+l,n)  =  A(k,n)A(k7n-2)  -  A(k,n-1)2   (28) 
A(k,n)+A(k,n-2)-2A(k,n-l) 

n=2k+3,2k+4, .... 
To  illustrate  the  power  of  this  computational  device  con- 
sider the  application  of  equation  (28)  and  its  iterates  to 
the  first  ten  terms  of  the  alternating  series 

Ln(2)  =  1"|+-|-|+---  (29) 

The  results  are 
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n      A(0,n)  A(l,n)  A(2,n)  A(3,n)  A(4,n) 


1  1.  oooccooooo  0.0 

2  0.5GG0COGOOQ  0.0 


where  the  tenth  partial  sum,  A(0,10),  is  correct  to  only  one 
significant  figure.   However,  iterate  A(4,10)  is  correct  to 
eight  figures  since  Ln(2)  =  0.6931471806. 

In  this  thesis  research,  Shanks'  accelerator  has  been 
tested  extensively  in  conjunction  with  the  simple  parameter- 
ized contours  of  integration  discussed  in  Chapters  II  and 
III. 

Modification  of  the  basic  algorithm  of  Chapter  III  con- 
sisted of  the  elimination  from  the  main  program  of  the  ter- 
mination criterion  for  the  numerical  integration  given  by 
equation  (16),  and,  instead,  performing  the  integration  for 
a  finite  number  of  terms.   These  addends  of  equation  (16) 
were  stored  in  a  linear  array  and,  subsequently,  transformed, 
within  the  main  program  in  accordance  with  equation  (28). 

The  accelerator  did  not  enhance  the  rate  of  convergence 
of  the  numerical  inversion  process  for  any  of  the  cases 
investigated.   In  fact,  if  the  number  of  terms  of  A  ,  which 
are  transformed,  becomes  sufficiently  large  so  that  the 
summation  of  the  addend  contributions  of  equation  (16) 
approaches  the  correct  result,  the  transformation  will 
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decrease  the  accuracy  of  the  result  as  the  denominator  of 
equation  (28)  becomes  small. 

Appendix  A  contains  the  program  listing  for  the  simple 
parameterized  contour  of  integration  in  conjunction  with 
Shanks'  accelerator. 

Shanks  [19]  discusses  transformations  of  higher  order. 
These  transformations  may  be  worthy  of  investigation. 
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VII.   LOCATION  OF  FUNCTION  SINGULARITIES 

The  criteria  of  Chapter  II  for  the  distortion  of  the 
Bromwich  contour  include  the  requirement  that  all  of  the 
singularities  of  a  function  f(s)  be  to  the  left  of  the  con- 
tour.  Heretofore,  the  discussion  has  assumed  that  the  loca- 
tion of  these  singularities  was  known.   Indeed,  this  is  not, 
in  fact,  always  the  case.   In  this  chapter  we  investigate  a 
method  of  determining  the  location  of  the  singularities  of 
f(s)  when  such  information  is  either  not  known  or  is  diffi- 
cult to  obtain  analytically. 

To  see  how  one  can  proceed  without  knowing  where  the 
singularities  are,  consider  the  simple  case 


f(s)  =  ~2~^ — 2  (30) 

s   +  a 


for  the  particular  values  of  a  =  0.125  and  t  =  2tt  . 

We  use  a  simple  parabolic  contour  given  by  equation  (13) 

with  B=l  and  with  A,  which  we  denote  by  the  symbol  x   in  the 

'  J  o 

discussion  to  emphasize  its  actual  significance,  being  given 

various  values.   If  we  plot  the  result  of  the  inversion 

versus  the  starting  position  x  ,  the  result  is  as  shown  in 

Figure  3. 

For  values  of  x   less  than  zero,  the  result  is  0.00. 

o 

For  x   in  the  range  from  about  0.2  to  2.0,  the  result  is 
0.707107.   For  x   larger  than  about  2.0,  the  result  is  near 
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0.7  but  is  sensitive  to  the  precise  value  of  x  .   For  x   in 
the  range  from  0.0  to  0.2,  the  result  exhibits  a  "spike" 
and  a  very  significant  change  from  0.0  to  0.707107. 
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Figure  3.   F(t)  vs  contour  starting  position  for 

cos(ctt).   Value  of  inverse  of  s_ 

for  t  =  2tt  as  a  function 


of  starting  abscissa  x  . 
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We  can  explain  this  behavior  as  follows.   For  x   less 
c  o 

than  zero,  the  singularities  are  on  the  wrong  side  of  the 
contour  and  we  simply  get  the  wrong  answer.   For  x   in  the 
neighborhood  of  0.0  to  0.2,  the  contour  passes  so  close  to 
the  singularity  that  the  trapezoidal  integration  scheme  is 
not  accurate  and  experiences  difficulty  in  getting  an 
accurate  value  either  for  the  right  answer,  0.707107  or  for 
the  wrong  answer  0.0.   For  x   in  the  range  from  0.2  to  2.0 
the  singularities  are  on  the  correct  side  of  the  contour  and 
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we  get  the  correct  answer.   For  x   larger  than  2.0,  there  are 
numerical  difficulties,  probably  due  to  oscillation  or  to 
the  effects  of  the  exp(st)  factor  for  s  having  large  real 
part. 

We  have  done  the  same  thing  with  many  other  transforms 
and  have  found  the  same  behavior.   Our  experience  leads  to 
this  suggestion  for  cases  where  the  user  does  not  know  the 
location  of  the  singularities.   Perform  the  calculation  for 
a  number  of  different  values  of  x   and  plot  the  results  as 
in  Figure  3.   The  correct  result  is  the  ordinate  of  the 
rightmost  plateau,  before  oscillations  or  other  numerical 
difficulties  have  had  a  chance  to  introduce  inaccuracies. 
In  the  case  illustrated  in  Figure  3,  there  were  only  two 
plateaus,  one  corresponding  to  all  the  singularities  lying 
on  the  wrong  side  of  the  contour  and  one  corresponding  to 
all  the  singularities  lying  on  the  correct  side  of  the  con- 
tour.  However,  with  more  singularities,  one  might  expect 
to  have  several  plateaus,  only  the  rightmost  of  which  repre- 
sents the  correct  result. 
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VIII.   THE  HYPERBOLIC  TANGENT 

In  Chapter  II  it  was  stated  that  great  advantage  accrues 
from  distorting  the  Bromwich  contour  as  long  as  the  process 
of  distortion  does  not  cause  the  contour  to  cross  over  a 
singularity  of  f(s).   In  this  chapter  we  investigate  a  typi- 
cal case  in  which  it  is  impossible  to  satisfy  this  condition 

The  function 


f  (s)  =  |  tanh(^) 


(31) 


has  an  infinite  sequence  of  simple  poles  on  the  imaginary 

IT 

axis  at  points  s  =  ±(2n+l)—  .   By  other  methods  it  may  be 
shown  that  this  function  is  the  transform  of  the  square  wave 
function  shown  in  Figure  4.   In  this  chapter  we  investigate 
the  extent  to  which  our  present  methodology  is  capable  of 
obtaining  this  square  wave. 
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Figure  4.   Square  wave  function. 

The  convergence  of  the  numerical  integration  has  been 
assured  by  bending  the  contour  to  the  left  in  all  other  cases 
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reported  in  this  thesis.   Accordingly,  we  do  so  in  the  pres- 
ent instance.   This  implies  that  we  must  pass  between  two 
poles  as  the  contour  bends  to  the  left  and  all  of  the  in- 
finite number  of  poles  above  the  point  where  the  contour 
crosses  the  imaginary  axis  thus  fall  on  the  wrong  side  of 
the  contour  and  are  not  included  on  the  left.   Thus,  in 
effect,  our  procedure  is  finding  the  inverse  Laplace  trans- 
form of  a  substitute  function  which  has  only  a  finite  number 
of  poles  coinciding  in  location  with  those  of  the  given 
function  f(s)  and  having  the  same  residues.   There  is  some 
reason  to  hope  that  the  inverse  of  the  substitute  function 
will  be  essentially  similar  to  that  of  the  given  function, 
at  least  for  sufficiently  small  values  of  t. 

Indeed  we  do  find  this  to  be  the  case.   SUBROUTINE  CURVE 
of  the  algorithm  of  Chapter  II  was  modified  to  formulate 
a  simple  parameterized  curve,  which  utilized  the  Bromwich 
contour  for  numerical  integration  in  its  initial  segment, 
starting  from  a  position  s=y  on  the  positive  real  axis  and 
continuing  until  the  elevation  of  a  finite  number  of  poles 
is  reached.   The  second  segment  of  the  simple  parameterized 
curve  which  was  used  consists  of  the  parabola  of  equation 
(18).    This  parabolic  arc  allowed  the  integration  contour 
to  cross  the  imaginary  axis  between  adjacent  poles  of  the 
function  f(s).   Such  a  contour  is  indicated  in  Figure  3. 

The  termination  criterion  which  was  used  in  this  case 
is  identical  to  that  described  in  Chapter  III.  All  other 
program  input  values  remained  unaltered  from  the  form  of  their 
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Figure  5.   Contour  of  integration  for  the  hyperbolic 
tangent..  Special  contour  for  inversion  of 
f(s)  =  ^  tanh(cts/2).   It  starts  at  x  , 
rises  vertically,  and  then  follows  a  para- 
bolic path  between  two  poles. 


previous  presentation  with  the  exception  that  the  integer 
JJ JJJ ,  which  selected  the  number  of  function  poles  to  be 
encompassed  by  the  contour  was  added.   The  program  listing 
appears  in  Appendix  A. 

The  approximate  solutions  F(t)  obtained  by  numerical 
inversion  of  the  substitute  function  for  the  f(s)  given 
by  equation  (31)  are  shown  graphically  below  in  Figures  6-9. 

Figure  6  is  a  reproduction  of  five  periods  of  the  square 
wave  of  Figure  4  with  a  =  10.0  and  the  first  100  poles  of 
f(s)  above  the  real  axis  to  the  left  of  the  contour  of  in- 
tegration.  Figure  7  is  a  similar  picture  of  the  same  F(t) 
with  a  =  5.0  and  the  first  150  poles  of  f(s)  above  the  real 
axis  encompassed  by  the  contour,  with  only  one  period  of  the 
wave  form.   Figures  8  and  9  are  analogous  to  Figures  6  and  7, 
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Figure  6. 


Numerical  approximation  of  five  cycles  of  the 
square  wave  function  with  100  poles  to  the 
left  of  the  contour. 
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Figure  7.   Numerical  approximation  of  one  cycle  of  the 
square  wave  function  with  150  poles  to  the 
left  of  the  contour. 
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Figure  8.   Numerical  approximation  of  five  cycles  of 
the  square  wave  function  with  7  poles  to 
the  left  of  the  contour. 
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Figure  9.   Numerical  approximation  of  one  cycle  of  the 
square  wave  function  with  7  poles  to  the 
left  of  the  contour. 
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respectively;  however,  only  the  first  seven  poles  of  f(s) 
above  the  real  axis  were  included  on  the  correct  side  of  the 
contour  of  integration. 

The  reader  should  be  aware  that  the  increment  of  t  used 
in  the  above  series  of  graphs  differed  in  each  case.   The 
inaccuracies  shown  in  Figures  6-9  are  partially  due  to  the 
assumptions  inherent  within  the  employment  of  our  algorithm, 
and,  also,  partially  due  to  the  fact  that  the  plotter  simply 
connects  points  sequentially.   This  can  also  be  seen  in  the 
slope  of  the  graphs  at  t  =  a ,  2a,  etc.,  which  is  a  function 
of  the  particular  At  employed  and  not  an  inaccuracy  of  the 
algorithm. 

Efforts  to  accomplish  numerical  inversion  of  the  f(s) 
of  equation  (31)  using  the  steepest  descent  contour  were 
unsuccessful.   The  difficulty  in  using  this  contour  appeared 
to  be  related  to  an  adverse  influence  from  the  proximity  of 
the  function  poles  to  the  v  =  0  contours  near  the  imaginary 
axis.   The  algorithm  was  overtaxed  in  this  vicinity  and 
could  not  track  the  contour;  consequently,  further  investi- 
gation was  not  attempted. 

In  summary,  we  conclude  that  the  effect  of  substituting 
a  function  which  has  only  a  finite  number  of  poles  which  are 
coincident  in  location  with  the  first  n  poles  of  a  function 
f(s),  which  possesses  an  infinite  number  of  poles  spaced  incre- 
mentally  along  the  imaginary  axis,  produces  an  inverse 
approximation  to  the  inverse  of  the  given  function  f(s). 
The  accuracy  obtained  with  such  a  substitution  is  increased 
as  more  of  the  function  singularities  are  encompassed  by  the 

distorted  contour. 
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IX.   FACTORS  WHICH  INFLUENCE  ACCURACY 

Heretofore,  our  discussion  has  included  the  treatment  of 
contour  integration,  utilizing  a  simple  parameterized  curve 
for  that  purpose,  from  two  perspectives.   The  first  of  these 
methods,  discussed  in  Chapter  III,  incorporated  a  termina- 
tion criterion  which  required  a  prescribed  number,  N,  of 
successive  integrand  contributions  each  to  have  smaller  mag- 
nitude than  a  specified  epsilon.   The  second  of  these  methods, 
discussed  in  Chapter  VI,  performed  numerical  integration  for 
a  finite  number  of  points  spaced  along  the  contour  with  sub- 
sequent employment  of  Shanks'  accelerator  in  an  endeavor  to 
enhance  convergence. 

The  first  of  these  methods,  we  have  found,  was  highly 
successful  when  employed  with  the  proper  point  spacing,  and 
a  very  tight  termination  criterion,  namely,  a  small  epsilon 
and  large  N.   The  second  method  has  not,  in  our  investiga- 
tions, been  found  to  be  successful,  and  need  not  be  discussed 
further.   The  issue  then  becomes  one  of  selecting  the  proper 
combination  of  step-size,  epsilon,  and  N  to  provide  an 
accurate  and  economical  result. 

In  order  to  make  such  an  appropriate  selection  of  these 
parameters,  it  is  necessary  to  examine  the  factors  which  in- 
fluence accuracy.   These  factors  include:  the  step-size  along 
the  contour,  the  termination  criteria,  and  the  oscillation. 
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Consider  a  simple  case,  which  we  have  chosen,  along  a 
path  for  which  the  oscillation  is  insignificant.   The  accu- 
racy in  such  a  case  is  influenced,  if  the  termination 
criteria  are  sufficiently  stringent,  only  by  the  step-size 
along  the  contour  of  integration.   Thus  dealing  with  the 
transform 

f(s)  =  -%-$ — ^  (32) 

s   +  a 

whose  inverse  is  the  cosine  function,  yields  a  plot  of  the 
logarithm  of  the  percent  error  of  F(t)  versus  the  step-size, 
A,  as  shown  in  Figure  10.   In  this  case,  the  contour  of  in- 
tegration was  the  simple  parameterized  curve  of  equation 
(18)  with  A  =  0.3,  B  =  1.0,  a  =  0.125  and  t  =  2tt. 

The  sequences  of  points  in  this  plot  are  all  obtained  by 
maintaining  N  constant  (N=3)  and  plotting  log-iQ  percent  error 
versus  A  with  epsilon  as  a  parameter. 

Four  such  parameteric  plots  appear  in  the  figure  below 
corresponding  to  the  values  of  epsilon  (EPS)  as  shown  in 
the  legend. 

Clearly  the  accuracy  of  the  results  may  be  seen  to  behave, 
in  a  general  sense,  in  the  manner  shown  in  Figure  11,  where 
EPS1  <  EPS2  <  EPS3  <  EPS4. 

When  the  termination  criteria  are  tight,  namely,  when  EPS 
is  small  and  N  is ,  at  least,  greater  than  one,  the  percent 
error  is  a  narrowly  banded,  roughly  linear  function  of  the 
integration  step  size,  A,  over  a  range  from  about  A  =  0.02 
to  A  =  0.2.   This  corresponds  to  EPS-,  in  Figure  11. 
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Figure  10.   Logarithm  of  percent  error  versus  step  size 
for  the  cosine  function. 


If  the  termination  criterion  EPS  is  looser,  as  evidenced 
in  other  cases,  shown  schematically  in  Figure  11,  the  accu- 
racy of  the  result  cannot  be  improved  indefinitely  by  decreas- 
ing, A.   Rather,  there  is  a  point,  which  differs  in  each  case, 
beyond  which  the  accuracy  is  decreased  as  the  step  size  is 
decreased.   The  behavior  is  complicated  by  the  fact  that  the 
effects  of  changes  in  step-size  and  changes  in  termination 
criteria  are  not  themselves  disjoint.   Making  the  step-size 
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Schematic  representation  of  accuracy  factors 
Behavior  of  the  logarithm  of  percent  error 
versus  step  size  A  as  a  function  of  the  ter- 
mination criterion  EPS. 


smaller  results  in  meeting  the  termination  criteria  at  an 
earlier  point  along  the  path  of  integration  and  thus  may 
actually  reduce  overall  accuracy. 

One  way  of  dealing  with  this  matter,  of  course,  is  to 
revise  the  termination  criteria  so  as  to  remove  this  inter- 
dependence.  We  simply  have  not  experimented  with  using 
alternate  forms  for  the  termination  algorithm.   Instead,  we 
suggest  the  following  viewpoint  to  the  user  who  wishes  to 
assure  that  he  obtains  results  with  a  specified  accuracy. 

He  should  first  of  all  employ  a  very  strict  termination 
criterion  and  use  a  succession  of  rather  small  values  of 
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step  size,  and  should  also  vary  the  starting  point  on  the 
contour.   In  this  way  he  should  be  able  to  establish  an 
evaluation  having  even  greater  accuracy  than  he  requires. 
Then,  selecting  what  seems  to  be  a  good  starting  point  and 
keeping  to  it,  he  should  increase  the  step  size  until  the 
result  is  no  more  accurate  than  he  requires.   He  should 
observe  how  much  saving  in  computer  time  (or  in  the  number 
of  calls  to  VALUE)  this  affords  and  should  not  increase  step 
size  beyond  a  reasonable  point.   Then  he  should  loosen  the 
termination  criterion  progressively,  measuring  the  saving 
in  computer  time  versus  the  possible  loss  of  accuracy, 
stopping  short  of  the  point  where  he  cannot  rely  on  obtain- 
ing the  required  accuracy. 

In  many  cases  where  numerical  integration  is  employed 
to  obtain  a  result  whose  correct  value  is  Y,  the  numerically 
produced  approximation  y  approximately  satisfies  a  relation 

y  =  Y  +  aAm  (33) 

where  a  and  m  are  unknown  constants.   The  exponent  m  may 
frequently  be  established,  once  and  for  all,  for  the  type 
of  integration  being  used,  and,  for  individual  integrands, 
the  correct  value,  Y,  from  two  evaluations,  y,  obtained  with 
A  =  A,  and  y2  obtained  with  A  =  A„,  by  using  the  extrapola- 
tion formula 


Y  "  17%  tv2  -  ViJ  (34> 
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We  have  investigated  this  method  of  obtaining  improved  accu- 
racy for  our  inversions  and  have  found  that  with  the  inter- 
dependence between  termination  criteria  and  step  size,  the 
numbers  a  and  m  were  both  rather  large  and  difficult  to 
establish  -  in  other  words,  the  extrapolation  was  not  success- 
ful in  producing  improved  values.   Our  conclusion  was  that 
it  was  more  efficient  from  a  computational  point  of  view 
simply  to  use  a  sufficiently  small  step  size  and  an  appro- 
priately matched  set  of  termination  criteria  so  as  to  be 
able  to  obtain  an  accurate  answer  without  interpolation. 

However,  a  suggestion  for  further  study  of  this  matter 
is  made  at  the  end  of  the  next  chapter. 
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X.   RESULTS,  CONCLUSIONS,  AND  RECOMMENDATIONS 

This  thesis  has  demonstrated  that  numerical  integration 
along  a  suitable  contour  in  the  complex  plane,  thus  implement- 
ing the  complex  inversion  integral  formula,  is  an  effective 
way  of  obtaining  numerical  values  for  the  inverse  Laplace 
transformation  of  a  given  function. 

The  success  of  the  procedure ,  i.e.,  its  efficiency,  avoid- 
ance of  inaccuracy  due  to  oscillation,  and  its  termination 
after  including  only  a  reasonable  number  of  addends  in  the 
numerical  sum,  is  related  to  the  fact  that  the  contours  chosen 
bend  to  the  left  as  they  rise,  thus  taking  advantage  of  the 
exp(st)  factor  in  the  integrand  of  the  inversion  formula. 

Although  there  are  some  theoretical  advantages  to  what 
we  have  called  the  "steepest  descent"  contour  which  follows 
a  v=0  path,  the  algorithm  which  permits  following  such  a  path 
requires  repeated  evaluations  of  the  integrand,  many  of  which 
are  not  actually  used  in  forming  the  summands  for  the  numeri- 
cal summation.   We  have  found  by  investigating  almost  100 
cases  for  which  known  inverses  are  available  in  analytic  form, 
that  it  is  more  efficient  to  use  a  simpler  contour  than  the 
steepest  descent  path.   By  using  what  we  have  called  a  simple 
parameterized  curve,  we  have  devised  a  procedure  which  makes 
use  of  every  evaluation  of  the  integrand.   We  have  found  that 
the  penalties  of  not  following  a  steepest  descent  curve, 
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namely  an  increase  in  oscillation  of  the  integrand  along  the 
path,  may  be  made  quite  tolerable  by  using  any  of  a  number 
of  simple  curves.   In  particular,  most  of  our  successful  work 
has  been  with  a  simple  parabola.   Our  algorithm,  however, 
includes  the  output  of  information  which  can  serve  to  alert 
the  user  to  the  possibility  of  inaccuracy  due  to  oscillation 
and  thus  suggest  to  him  that  he  might  select  a  different  form 
for  the  path. 

We  have  even  been  successful  in  obtaining  the  inverse  of 
a  function  which  has  an  infinity  of  poles  spaced  along  the 
imaginary    axis.   In  doing  so  we  violated  the  principle 
that  all  poles  of  the  transform  must  be  to  the  left  of  the 
contour.   Nevertheless  the  results  are  quite  satisfactory. 

Our  success  and  accuracy  have  been  so  gratifying  that 
we  venture  to  suggest  that  our  method  may  prove  to  be  an 
efficient  alternate  method  for  the  evalution  of  exotic  func- 
tions for  which  other  methods  are  slowly  convergent  or  in- 
volve series  the  coefficients  in  which  are  difficult  to 
obtain.   For  example,  case  92  in  Appendix  B  shows  the  suc- 
cessful evaluation  of  the  rather  uncommon  Struve  function. 

The  original  reason  for  attempting  to  employ  numerical 
integration  in  the  complex  plane  as  a  means  of  inverting  a 
Laplace  transformation  was  that  alternate  methods  employed 
by  Hiep  and  Zargary  in  conjugated  heat  transfer  problems 
were  simply  not  accurate  enough.   They  led  to  physically 
impossible  results  with  some  temperatures  in  the  media  below 
sink  temperature  and  others  above  source  temperature.   In 
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this  thesis  we  have  not  restudied  these  heat  transfer  prob- 
lems and  recommend  that  this  be  done  using  our  methods  of 
inversion.   There  is  some  reason  to  believe  that  the  diffi- 
culty encountered  by  Hiep  and  Zargary  will  not  be  encountered 
if  our  method  is  used,  but  we  are  not  prepared  to  substan- 
tiate this  claim  at  this  time. 

Also  it  might  be  of  use  if  the  effects  of  varying  the 
available  parameters  (shape  of  curve,  spacing  of  points, 
starting  point  of  curves,  termination  criteria,  etc.)  were 
to  be  investigated  further  so  that  a  user  would  be  better 
prepared  to  deal  with  indications  of  inaccurate  inversion 
or  inordinate  requirements  for  computer  time.   Our  own 
numerical  experiments  were  invariably  so  happily  successful 
that  we  have  not  encountered  need  for  such  information. 

For  functions  f(s)  for  which  it  may  be  difficult  to 
locate  the  singularities,  we  have  shown  that  varying  the 
starting  point  of  a  simple  parameterized  contour  permits 
selecting  an  optimum  contour  in  the  sense  that  one  can  be 
assured  that  all  singularities  are  to  the  left  of  the 
curve  and  also  that  the  curve  does  not  reach  far  enough  to 
the  right  to  impose  numerical  difficulties  with  large  posi- 
tive exponentials. 

At  the  end  of  the  preceding  chapter  we  indicated  that 
we  did  not  find  it  profitable  to  employ  extrapolation  as  a 
means  of  obtaining  improved  accuracy.   However,  this  is 
probably  worth  looking  at  again,  and  our  suggestions  for 
doing  so  are  as  follows.   First  employ  a  termination  criterion 
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which  is  disjoint  from  the  integration  step  size  in  the 
sense  that  the  location,  along  the  path  of  integration,  of 
the  point  where  termination  occurs,  is  independent  of  step 
size.   One  way  of  doing  this  is  to  employ  an  epsilon  in  the 
termination  criterion  which  is  itself  the  result  of  multiply- 
ing a  fixed,  input  epsilon  by  the  step  size.   Then  the 
extrapolation  has  a  chance  of  being  successful.   So  as  to 
maintain  optimum  computational  efficiency,  one  should  use 
what  is  called  Richardson  extrapolation  in  which  the  evalua- 
tion points  for  the  larger  step  size  are  themselves  included 
among  those  for  the  smaller  step  size. 
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APPENDIX  A 
LISTING  OF  COMPUTER  PROGRAMS  AND  SUBROUTINES 

The  computer  programs  and  principal  subroutines  that  we 
have  prepared  and  tested  during  our  investigation  are  all 
listed  within  this  appendix.   Section  A-l  contains  the  user's 
instructions  glossary,  flowchart,  and  program  listing  for  the 
simple  parameterized  curve.   Section  A-2  is  a  similar  treat- 
ment for  the  steepest  descent  contour.   Section  A-3  contains 
the  user's  instructions  and  program  listing  for  the  simple 
parameterized  curve  procedure  adapted  for  use  with  Shanks' 
accelerator.   Section  A-4  contains  the  same  information  for 
adaptation  of  the  simple  parameterized  curve  for  use  with 
the  hyperbolic  tangent  function  described  in  Chapter  VIII. 

It  is  hoped  that  sufficient  detail  has  been  provided  to 
enable  a  user  to  adapt  one  of  these  programs  to  his  purpose 
in  an  efficient  manner.   However,  if  additional  insight  is 
required,  it  may  be  necessary  to  refer  to  the  chapter  of 
this  thesis  in  which  the  algorithm  is  developed.   This  is 
particularly  true  in  the  cases  of  Sections  A-3  and  A-4;  the 
development  in  these  sections  has  not  been  repeated  where  it 
is  equivalent  to  that  of  Section  A-l. 

Functions  DREAL  and  DIMAG  and  SUBROUTINES  CPOWER  and 
CLOG,  if  they  are  required,  are  listed  within  Chapter  III, 
as  is  a  typical  example  of  SUBROUTINE  VALUE.   Additionally, 
all  functions  f(s)  which  have  been  investigated  during  this 
research  have  their  applicable  SUBROUTINE  VALUE  listed  within 
Appendix  B. 
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SECTION  A-l:   THE  SIMPLE  PARAMETERIZED  CURVE 

USER  INSTRUCTIONS 


This  is  the  FORTRAN  IV  program  for  the  numerical  inver- 
sion of  the  Laplace  transformation  of  a  function  f(s),  which 
performs  contour  integration  along  a  distorted  Bromwich  con- 
tour in  the  form  of  a  simple  parameterized  curve.   The  fol- 
lowing instructions  are  provided  to  assist  the  user  in 
familiarizing  himself  with  the  program  so  that  he  can  adapt 
it  to  his  particular  requirements. 

1.  SUBROUTINE  CURVE  is  a  user  supplied  group  of  instructions 
which  formulates  the  numerical  integration  contour.   The 
program  listing  presently  contains  a  simple  parabola,  which 
may  be  utilized  directly  if  desired. 

2.  SUBROUTINE  VALUE  is  a  user  supplied  subroutine  which 
calculates  the  real  and  imaginary  components  for  g(s)  = 
exp(st)*f(s)  =  u(s)  +  iv(s).   For  examples  of  preparation 
of  this  subroutine  the  user  is  referred  to  the  test  cases 
of  Appendix  B. 

3.  The  mandatory  input  variables  to  the  main  program  are 
as  follows: 

(a)  AA  -  REAL*8  starting  position  (S=AA)  on  the  real 
axis  of  the  contour  of  integration.   (Also  denoted 
by  A  in  equation  (18)  and  x   in  Chapter  VII.) 

(b)  AL  -  REAL*8  value  of  any  constant  associated  with 
the  function  f (s).   (Note:   If  the  function  f(s) 
involves  more  than  one  constant,  the  calling  state- 
ments to  SUBROUTINE  VALUE  must  be  modified  accord- 
ingly. 
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(c)  DELP-REAL*8  increment  of  the  p-parameter  leading  to 
spacing  of  points  along  the  contour  of  integration. 

(d)  T  -  REAL*8  time  for  which  F(t)  is  desired. 

4.  Input  parameters  which  may  be  altered  by  the  user,  at 
his  discretion,  are  as  follows: 

(a)  EPS  -  REAL*8  tolerance  parameter  for  the  termina- 
tion of  numerical  integration.   The  program  listing 
presently  contains  a  value  of  l.D-ll. 

(b)  NUMBER  -  INTEGER*4  number  of  successive  addend  con- 
tributions less  than  the  specified  value  of  EPS 
for  which  numerical  integration  is  terminated. 

5.  The  output  variables  which  are  printed  as  output  by  pro- 
gram are  defined  below: 

(a)  AA     starting  position  on  the  real  axis 

(b)  SUM    the  value  of  F(t)  obtained  by  numerical 

inversion  of  f(s).   (cf.  equation  (16)) 

(c)  SUMA    the  absolute  value  of  the  summation  of  the 

the  addend  contributions  of  the  numerical 
integration,  (cf.  equation  (14)) 

(d)  X      the  final  value  of  x  on  the  contour  at  which 

integration  was  terminated 

(e)  Y      the  final  value  of  y  on  the  contour  at  which 

integration  was  terminated 

(f)  LOSC   the  number  of  sign  changes  between  successive 

addends  encountered  during  numerical  integra- 
tion 

(g)  MMMMM   the  total  number  of  calls  made  to  SUBROUTINE 

VALUE  during  the  numerical  integration  process 


6.  The  output  of  these  variables  is  as  follows: 

AA,  SUM,  SUMA,  X,  Y,  LOSC ,  MMMMM 

7.  The  time  (seconds)  required  to  perform  the  numerical  inte- 
gration is  printed  on  a  line  immediately  preceding  the  items 
listed  above. 

8.  As  previously  described,  the  program  will  perform  numeri- 
cal integration  along  one  simple  parameterized  curve.  If  the 
user  desires,  the  program  may  be  modified  to  "sweep"  the  con- 
tour of  integration  over  a  range  of  starting  values,  AA,  (cf. 
Chapter  VII).  This  may  be  quickly  accomplished  in  the  follow- 
ing manner. 

(a)  Replace  the  value  of  the  integer  KSTART  from  1  to 
the  number  of  contour  integrations  desired. 

(b)  Replace  the  REAL*8  value  of  AA  with  the  following 
statement 

A A  =  XSTART+AK*XINCRM 

where 

XSTART  =  the  smallest  starting  value  of  the 
contour  of  integration  the  user 
desires 

XINCRM  =  the  Ax  between  starting  positions 

of  successive  contours  of  integration 
must  also  be  supplied. 

Thus,  the  program  will  perform  numerical  inversion  along 
KSTART  contours  of  integration  which  differ  only  in  their 
starting  positions  along  the  x  axis.   The  user  may  then 
observe  increases  or  jumps  in  the  value  of  the  output  variable 
SUM  between  successive  integrations,  enabling  him  to  ascertain 
the  locations  of  function  poles.   The  utility  of  this  proce- 
dure is  discussed  in  depth  in  Chapter  VII. 
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9.   If  SUBROUTINES  CPOWER  and/or  CLOG  are  called  by  SUBROU- 
TINE VALUE  in  conjunction  with  use  by  this  program,  NCALL 
and/or  MCALL  must  be  dimensioned  and  initialized  to  zero 
within  the  main  program  and  added  to  the  calling  statements 
to  SUBROUTINE  VALUE.   Additionally,  any  powers,  other  than 
integer  values,  must  be  entered  within  the  main  program  and 
passed  within  the  calling  statements  to  SUBROUTINE  VALUE. 
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VARIABLE 
AA 

ADD 

ADDA 

ADDDLD 

AL 
DELP 

EPS 

KOUNT 

KSTART 


LOSC 


MMMMM 


NN 


NUMBER 


PROD 


GLOSSARY  OF  VARIABLE  NAMES 

DESCRIPTION 

starting  position  on  the  real  axis 
of  the  simple  parameterized  curve 

computed  value  of  the  addend  of  the 
numerical  integration 

absolute  value  of  the  addend  of  the 
numerical  integration 

saved  value  of  the  previous  addend 
contribution 

constant  associated  with  f(s) 

increment  of  point  spacing  along  the 
simple  parameterized  curve 

tolerance  parameter  for  termination 
of  the  numerical  integration  process 

integer  counter  utilized  to  terminate 
the  numerical  integration  process 

integer  counter  for  the  number  of  times 
the  numerical  integration  is  desired 
to  be  performed,  using  different  start- 
ing positions 

integer  counter  utilized  to  record  the 
number  of  sign  changes  between  succes- 
sive addends 

integer  counter  to  record  the  total 
number  of  calls  made  to  SUBROUTINE  VALUE 
during  the  numerical  integration 

integer  counter  to  record  the  total 
number  of  addend  contributions 

integer  input  of  the  desired  number  of 
successive  addend  contributions  less 
than  a  specified  epsilon  for  termination 

parameter  of  the  simple  parameterized 
curve 

flag  utilized  to  detect  oscillation  of 
the  integrand 


LEGEND 
2 


2 

2 

4 
1 
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SUM         total  addend  contribution  for  the  3 

numerical  integration  process 

SUMA        absolute  value  of  the  total  addend         3 
contribution  for  the  numerical 
integration  process 

T  time  for  which  F(t)  is  to  be  evaluated     2 

TEST        flag  for  termination  of  the  numerical      1 
integration  process 

U  real  part  of  the  function  g(s)  1 

UO  previous  value  of  U  1 

V  imaginary  part  of  the  function  g(s)  1 
VO  previous  value  of  V  1 
X  real  part  of  s  1 
XO  previous  value  of  X  1 

Y  imaginary  part  of  s  1 
YO  previous  value  of  Y  1 

LEGEND 

1.  No  action  required  by  the  user 

2.  Mandatory  user  input 

3.  Appears  as  program  output 

4.  May  be  altered  by  the  user  at  his  discretion 
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ESTABLISH  THE 

SYSTEM  LIMITS 

AND  PARAMETERS 

INITIALIZE 
SET  THE  TIMER 


CALL 
CURVE  AND  VALUE 
OBTAIN 
INITIAL 

U  AND  V 


DO  400 

INCREMENT  THE 

STARTING  POINT 

COMMENCE 

NUMERICAL 

INTEGRATION 
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i 

200 
CALL 

INCREP 
CURVE 
VALUE 


COMPUTE  THE 

INTEGRAND 

AND  ITS 

ABSOLUTE 

VALUE 


SUM  THE 
CONTRIBUTIONS 

OF  THE 
INTEGRAND 

AND  ITS 
ABSOLUTE  VALUE 


IF  NN 


EQ.  1 


ESTABLISH 

THE  SIGN  OF 

THE  INTEGRAND 


<- 
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COMPUTE  THE 

PRODUCT  OF  THE 

OLD  AND  NEW 

INTEGRANDS 


LOSC 


IF 
PRODUCT 


^T.O.D+0 


M 


THE  COUNTER 

FOR  THE 

NUMBER  OF 

SIGN  CHANGES 

IS  INCREMENTED 


KOUNT 


IF 
TEST 


LT.EPS 


THE  COUNTER 
FOR  INTEGRAND 

MAGNITUDE 
IS  INCREMENTED 


IF 
KOUNT 


LT. NUMBER 


■>GO  TO  200 
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STOP  THE 
TIMER 


WRITE 

THE 
OUTPUT 


A. 


400 

CONTINUE 
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C               THIS    IS    THE    PROGRAM    TC    PERFORM    NUMERICAL    INVERSION    CF    TFE    LAPLACE  10 

C              TRANSFORMATION    OF    A    FUNCTION    USING    THE     SIMPLE     PARAMETERIZED    CURVE  20 

IMPLICIT    REAL*3    I A-H, C-Z)                                                                           i    C^P\\-  30 

CALL    ERRSET     (  ^07,  256,-  1,  1, 1 ,  2  OS)                                                        LL'CQ  .   j  40 

NUMBER    =    5 *-  £.                                      50 

KSTART     =    1 »  60 

TWOPI    =    4.0+0*0ARSIN(  1.0+0)  ^                                   70 

AL   =    1.250-1 >~  9                                     80 

C  90 

C                   4X4*444»44ti44444444444444U<i444444444444«4444444444444*4*44*4.*4ti44d  100 

C             #  *              110 

C               *       THE    OUTER    <_00P    MAY    8E    UTILIZED    TC    SWEEP    THE    CCNTCUR    STARTING       *  120 

C              ((POSITION    OVER    A    OESIREO    RANGE    AS    DESCRIBED    IN     ThE     USER    INST.  *               130 

C              *  *               140 

C  150 

00    400     KKK=1, KSTART  160 

ADO    =     O.D+0  170 

ADDA    =     C.D+O  180 

SUM   =    O.C+0  190 

SUMA    =     C.D+O  200 

LOSC    =    0  210 

MfrMMM    =    0  22G 

AK   =    KKK  230 

AA   =    3.0-1 240 

T   =    TWOPI    >-  h                                  250 

P    =    O.D+0  -                                 260 

NN   =    0  270 

KCUNT    =0  .                                  280 

EPS   =     1.D-1L _  d,                                 290 

DELP    =     1.0-] ^  300 

CALL    CURVE     (AA,P,XG,YC)  _                                   210 

CALL    VALUE     (  XO  ,  YD  ,  T  ,A  L,  UO,  VO,  MMMMMJ  320 

TIME    =    O.D+0  330 

CALL    SETIME  ?40 

C  -30 

C                 $i$$5i$5$$<!S$$I5S5SSS$$5$SS$S$iSSSS$SSi5$iS$<SSii,$$iS$$S$iiSiiSiS  36  0 

C              $  $               370 

C      $   THE  FOLLOWING  LOOP  PERFORMS  THE  NUMERICAL  INTEGRATION  $      380 

200  NN  =  NN  +  1  390 

CALL    INCFEP     (P.GELPi  400 

CALL     CURVE    lAAtP.X,*}  410 

CALL    VALUE     (X  ,Y  ,  f  »AL  tU  *V  ,MMMMM  J  420 

ADD    =     <  I  V+VO)*<  X-XO)*<  U*UO)*(Y-YO)) /TWOPI  430 

IF     (MN.EO.l)     ADOOLD=AOD  440 

PROD    =     ACL.*ADDDLD  450 

IF    I  PROD. LT. C.D+O)     LGSC=LOSC+l  460 

ACDA    =     CABS(AOD)  470 

SUM    =    SUM+ADO  480 

SUMA    =     SUMA+ADOA  490 

C  500 

C               *  *                520 

C              *      THIS     IS    THE    TERMINATION    CRITERION    SEGMENT    OF     THE    ALGORITHM            *  530 

C             *  *              540 

TEST    =     DABS (ADD)  550 

IF    (TEST.LT.EPS)    K0UNT*KQUNT+1  560 

IF    (TEST. GE. EPS)     KOUNT-0  570 

CALL    GETI^E    (  IETJ  580 

EL   =    DFL0AT1IET  J*2.60-5  590 

TIME    =     TI^E+El  600 

CALL    SETIME  610 

IF    (KOUNT.EO. NUMBER!     GO    TO    300  62  0 

C               *  «               630 

C  6  50 

XO   =     X  660 

YO   =    Y  670 

UC    =    U  680 

VC    =    V  690 

AODOLD    =    AOD  700 

GO  TO  200  710 

C      $  $      720 

C       iJSSSSJ<$$$S$SiSiiJ*SJS$iS5SS$iiSJ$$$ii$iiS$SiSI<SJiJSiJ$iJ$$ii$4S  730 

C  740 
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300    CONTINUE  750 

WRITE     (to, 500)    NN.TIME  760 

WRITE     (to, 600)    AA» SUM,SUMA,X,Y,LOSC,MMMMM  770 

400    CONTINUE  780 

C  790 

C               *                                                                                                                                                                           *  800 

C                 4n*4*4*titi4),ii4**44<t***n*HH44****U#%Hn»**1in*nn****tiit*«444*4**4*ii*m  810 

C  820 

STOP                                      .  ■  830 

C  840 

C  850 

500   FCRMAT    (IX, 'TOTAL    ELAPSEC    TIME    AFTER    • ,  15 ,  •     ITERATIONS     IS     ',E20.5)  86C 

600    FORMAT     (  IX,  IP  5E  20  .5,  2  I  8  , /)  870 

END  880 

SUBROUTINE    CURVE    (A,P,X,Y)  10 

IMPLICIT    REAL*8    (A-H,0-Z)  20 

X    *    A-( P*P)  30 

Y   =    P  40 

RETURN  50 

END  60 

SUBROUTINE     INCREP     (P,0£LF)  10 

IMPLICIT    REAL*d    (A-M.O-Z)  20 

P    =    P*OELP  30 

RETURN  40 

END  50 
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SECTION  A-2:   THE  STEEPEST  DESCENT  CONTOUR 

USER  INSTRUCTIONS 


This  is  the  FORTRAN  IV  program  for  the  numerical  inver- 
sion of  the  Laplace  transformation  of  a  function  f(s),  which 
performs  contour  integration  along  the  steepest  descent  path 
in  the  complex  plane.   The  following  instructions  are  pro- 
vided to  assist  the  user  in  familiarizing  himself  with  the 
program  so  that  he  can  adapt  it  to  his  particular  require- 
ments. 

1.  SUBROUTINE  XMARCH  is  a  subroutine  which  tracks  a  speci- 
fied v=0  contour  (steepest  descent)  in  the  complex  plane, 
once  this  contour  has  been  intersected  for  the  first  time. 
As  shown  in  the  program  listing,  the  following  two  equations 

X  =  XHOLD  -  l.D-l*DCOS( THETA) 

Y  =  YHOLD  +  l.D-l*DS IN (THETA) 
march  the  numerical  integrations  in  the  favorable  direction 
along  the  contour.   The  equations  are  parametric  in  THETA, 
which  is  the  search  angle  provided  by  SUBROUTINE  ROTATE. 
The  l.D-1  is  a  polar  radius  which  may  be  altered  to  a  larger 
or  smaller  increment  by  the  user,  or  accepted  as  it  appears. 

2.  SUBROUTINE  ROTATE  increments  the  angular  search  angle 
THETA  in  order  to  locate  a  new  point  on  the  v=0  contour.   The 
angle  is  presently  incremented  by  5.D-2  radians  with  each 
call  to  the  subroutine.   This  may  be  altered  as  required. 

3.  SUBROUTINE  VALUE  is  a  user  supplied  subroutine  which  cal- 
culates the  real  and  imaginary  components  for  g(s)  =  exp(st) 

x  f(s)  =  u(s)  +  iv(s).   For  examples  of  preparation  of  this 
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subroutine,  the  user  is  referred  to  Appendix  B. 

4.  The  mandatory  input  variables  to  the  main  program  are 
listed  as  follows: 

(a)  AA  -  REAL*8  starting  position  (s=AA)  on  the  real  axis 

(b)  AL  -  REAL*8  value  of  any  constant  associated  with 

the  function  f(s).   (Note:  If  the  function  f(s) 
involves  more  than  one  constant,  the  calling 
statements  to  SUBROUTINE  VALUE  must  be  modified 
accordingly. ) 

(c)  DELY  -  REAL*8  increment  of  the  point  spacing  along 

the  Bromwich  contour.  (Note:  to  alter  the  point 
spacing  along  the  v=0  contour,  the  user  is 
referred  to  1  above.) 

(d)  NREG  -  INTEGER*4  input  value  which  specifies  the 

first,  second,  third,  etc.  v=0  contour  above 
the  positive  real  axis,  which  is  to  be  followed 
in  the  numerical  integration  process. 

(e)  T  -  REAL*8  time  for  which  F(t)  is  to  be  evaluated. 

5.  Input  parameters  which  may  be  altered  by  the  user,  at 
his  discretion,  are  listed  below  as  follows: 

(a)  EPS  -  REAL*8  tolerance  parameter   for  the  termina- 

tion of  numerical  integration.   The  program 
listing  presently  contains  a  value  of  l.D-11. 

(b)  NUMBER  -  INTEGER*4  number  of  successive  addend 

contributions  less  than  the  specified  value  of 
EPS  for  which  numerical  integration  is  terminated 
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6.  The  output  variables  which  are  received  from  the  program 
are  defined  below: 

(a)  MMMMM  -  the  number  of  calls  made  to  SUBROUTINE  VALUE 

(b)  SUM   -  the  total  addend  contribution  for  the  numeri- 

cal integration 

(c)  SUMA  -the  absolute  value  of  the  total  integrand 

contribution  for  the  numerical  integration 

(d)  EL    -  the  computation  time  in  seconds  required  for 

the  process 

7.  A  sample  program  output  is  as  follows 

694  MMMMM    \  ,         .   , 

0.724724D+0  0.324972D+1     SUM, SUMA/  a~  ^na  OI  vertlcal 

0.707107D+0  0.326734D+1     SUM,SUMAi  S;Ctl°n  , 
4027  MMMMM    f  at  en?  J1 

0.71090D+1  EL       i  computation 
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VARIABLE 
AA 

ADD 

ADDA 

AL 
CHECK 

DELY 

EPS 
EXAM 


KOUNT 


KSIGN 


LEVEL 


MMMMM 


NREG 
NUMBER 


GLOSSARY  OF  VARIABLE  NAMES 

DESCRIPTION 

starting  position  on  the  real  axis  of 
the  Bromwich  contour  segment  of  the 
integration 

computed  value  of  the  addend  of  numeri- 
cal integration 

absolute  value  of  the  addend  of  numeri- 
cal integration 

constant  associated  with  f(s) 

flag  to  ascertain  when  a  v=0  contour 
crossing  has  occurred 

increment  of  point  spacing  along  the 
Bromwich  contour  segment  of  numerical 
integration 

tolerance  parameter  for  termination  of 
numerical  integration 

flag  to  ascertain  when  a  contour  crossing 
has  occurred 

integer  counter  used  internally  for  test- 
ing contour  regions 

integer  counter  utilized  to  terminate 
the  numerical  integration  process 

integer  counter  utilized  to  record  the 
number  o  f  v=0  contour  levels  crossed 

integer  flag  assigned  to  regions  of 
positive  and  negative  v 

integer  counter  for  the  number  of  calls 
made  to  SUBROUTINE  VALUE  during  each 
segment  of  the  numerical  integration 
process 

the  desired  v=0  contour  level  above  the 
positive  real  axis  which  is  to  be  tracked 

integer  input  for  the  number  of  succes- 
sive addend  contributions  less  than  the 
tolerance  parameter  required  for  termi- 
nation 


LEGEND 
2 


2 

1 
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SUM 
SUMA 

T 
TERM 

TEST 

THETA 

THETAO 

THNEW 

THONE 

THTWO 

U 
UHOLD 

V 

VHOLD 

VNEW 

VOLD 

X 

XHOLD 

Y 

YNEW 


total  addend  contribution  for  the  3 

numerical  integration  process 

absolute  value  of  the  total  addend         3 
contribution  over  the  contour  of 
integration 

time  for  which  F(t)  is  to  be  evaluated     2 

termination  parameter  for  the  absolute     1 
value  of  the  addend 

Flag  to  ascertain  when  a  contour  1 

crossing  has  occurred 

search  angle  used  in  the  contour  1 

tracking  algorithm 

previous  value  of  THETA  1 

predicted  angular  location  of  a  new        1 
point  on  the  contour  from  a  known 
previous  point  • 

coordinate  utilized  in  linear  inter-       1 
polation 

coordinate  utilized  in  linear  inter-       1 
polation 

real  part  of  g(s)  1 

previous  value  of  U  1 

imaginary  part  of  g(s)  1 

previous  value  of  V  1 

actual  value  of  the  imaginary  part  of      1 
g(s)  for  the  predicted  v=0  location 

previous  value  of  V  1 

real  part  of  s  1 

previous  value  of  X  1 

imaginary  part  of  s  1 

predicted  elevation  of  the  desired  v=0     1 
contour  level 
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LEGEND 

1.  No  action  required  by  the  user 

2.  Mandatory  user  input 

3.  Appears  as  program  output 

4.  May  be  altered  by  the  user  at  his  discretion 
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START 


ESTABLISH   THE 
PARAMETERS 
INITIALIZE 
SET   THE 
TIMER 


CALL 

CURVE  AND  VALUE 

OBTAIN  THE 
INITIAL  U  AND  V 


200 

INCREMENT  Y 

COMMENCE 

VERTICAL  MARCH 
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CALL 
CURVE  AND  VALUE 


COMPUTE  THE 

VALUE  OF  THE 

INTEGRAND  AND 

ITS  ABSOLUTE 

VALUE 


SUM  THE 

CONTRIBUTIONS 

OF  THE 

INTEGRAND 

AND  ITS 

ABSOLUTE  VALUE 


V 


SIGN 


IF 
J 


EQ.  2 


ESTABLISH 
THE  SIGN 
OF  V 


-<- 
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COMPUTE  THE 
PRODUCT  OF  THE 
OLD  AND  NEW  V 


K  SIGN 


IF 
PRODUCT 


LT.OD+0 


■> 


THE  COUNTER 

FOR  THE 

NUMBER  OF 

CONTOUR  LEVELS 

IS  INCREMENTS 


<■ 


300 
A- 


IF 
KSIGN 


LT . NREG 


JL 


CONDUCT  LINEAR 
INTERPOLATION 

LOCATE 

DESIRED  POINT 

ON  THE 

V=0  COUNTER 
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RETAIN  THE 
X,  Y,  U,  AND  V 
VALUES  OF 
THE  POINT 
ON  THE  COUNTER 


COMMENCE  THE 
TRACKING  ROUTINE 
TO  LOCATE  NEW 
POINTS  ON  THE 
COUNTER 


1100 
CALL 
X  MARCH 
ROTATE 

AND 
VALUE 


IF 
J 


EQ.  2 


ESTABLISH 

THE  SIGN 

OE  V 

SET  OLD  V  =  V 


■4 


Y 


74 


EXAM 
COMPUTE  THE 
PRODUCT  OF 
OLD  V  AND  V 


1100 


GT.O.D+0 


IF 
EXAM 


LT.O.D+O 


1200 
1300 


CALL 
VALUE 
OBTAIN  ACTUAL 
V  NEW  FOR  THE 
INTERPOLATED  9 


COMPUTE 
CHECK=V  NEW  *  SIGN 
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LT . EPS 


300 


REPLACE   THE 

COORDINATE 

PAIR  ASSOCIATED 

WITH  POSITIVE  V 

RE  INTERPOLATE 


IF 
DABS(V) 


GT.EPS 


1300 


GT.O.DH 


IF 
CHECK 


IT.O.D+0 


REPLACE   THE 

COORDINATE 

HPAIR  ASSOCIATED 

WITH  NEGATIVE  V 

RE INTERPOLATE 


1600 

COMPUTE  THE 

INTEGRAND 

AND  ITS 

ABSOLUTE  VALUE 


SUM  THE 

INTEGRAND 

AND  ITS 

ABSOLUTE  VALUE 
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1000 


COMPUTE 
TERM=DABS(ADD) 


IF 

TERM 


<r 


KOUNT 


.LT.EPS 


INCREMENT 
COUNTER 


IF 


LT.5 


KOUNT 

\          / 

\    EQ0 

5 

/ 

> 

f 

WRITE  OUTPUT 

STOP  THE  TIMER 

STOP 
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TRANSFORM    OF    A 


400 


5  00 


THIS    PROGRAM    COMPUTES    THE    INVERSE    LAPLACE 
FLNCTION    USING    Trit    STEEPEST    DESCENT    CONTOUR 
IMPLICIT    REAL*8     (A-H.Q-Z) 
CALL    SETIME 
CALL    ERRSET    (207,256,-1,1,1,2  091 

************* **********<**>)t*******$44*** 4  ********************£*£*,>: 

*  * 

*  THIS  SEGMENT  OF  THE  PROGRAM  NUMERICALLY  INTEGRATES  ALONG  """hE  * 
*eROM»VICH  CCNTOUR  UNTIL  IT  INTERSECTS  *  iTH  THE  UESIBEO  V»0  * 
•CONTOUR.       LINEAR    INTERPOLATION    1<    USEO    TO    REFINE    TFE    LCCATICN.       * 


ThOPI     =    4.0+0*QARSIN( 1  .C+0) 

T    =    TWOPI 

KKKKK    =    0 

AL    =    1.250-1 

J    =     1 
MMMM    =    0 
SLM    =    0.0+0 
SUM!,    =    0.0  +  0 

AA    =    3.C-1 

X    =    AA 

Y  =    O.D  +  0 

YFOLD    =    O.D+0 

OELY    =    2.50-3 

NREG    =     3  

KSIGN    = 


LE 


V^l 


o 


200 


3  00 


LEVEL  =  -I 
CALL  VALUE 
Y  =  Y+OELY 
CALL  VALUE 
J  =  J+l 
IF  (J.EC.2J 
IF    (J.GT.2) 


(X, Y, T,AL, LHCLD , VHOLD , f MPPM) 
(X,Y,T,AL,L,\.,MMMMM) 


SIGN    =    V/CA6S(V) 

TEST    =    V*SIGN 
IF    (TEST  .lE.O  .J+0  )    GO    TO    400 
ACO    =     < (U+UHOlOJ* (Y-YhCLCJ i/TWOPI 
SLM    =    SUM-t-ADO 
ACDA    =     CA8S  (ADO  J 
SUMA     =     SUMA+AOGA 

IF  (KKKKK.EC.U  WRITE  (6,1800)  ADO,AOOA 
UHOLD  =  U 
VHOLO  =  V 
YhOLD  =  Y 
GC  TO  2  00 
KSIGN  =  KSIGN+1 
LEVEL  =  -LEVEu 
IF  (KSIGN. cG.NREl,J 
SIGN  =  -SIGN 
GC   TO    300 


GO    TC    500 


YCNE     =     Y-OELY 

Y7*3    =     Y 

VCNE    =    VHOLD 

V7W0    =    V 
600    YNEW    =     (  VONE*(  YTWO-YONEJ /( VONE-VTViG) ) +YCNE 

CALL    VALUE     (X»Y>JEWtTt  ALfUNEW,  VNEW.MMMM.M) 

IF    (OAbS(VNEW).LT.l.G-UJ     GO    TO    9C0 

TEST    =     VNEW*SIGN 

IF    (TEST. GT  .0.0+0  J    GO    TO    700 

YTWO    =     YNEW 

V7W0    =     VNEW 

GO    TO    800 
700    YONE    =     YNEW 

VONfc    =     VNEW 
800    CCNTINUE 

GO    TO    6C0 
900    ACO   =     (  (UN  EV,*UHOLDJ*(  YNEW-YHOLO  )  ) /TUOPI 

ACOA    =     CA6S (AGO  J 

SLM    =     SLM+AOO 

SUMA    =     SUMA+AOOA 

IF    (KKKKK.EC.i)     WRITE     (6,1800J    AOCADQA 

WRI  TE     (  6,2100)     MMMMM 

WRITE     (6,2000)     SUM, SUMA 


10 

20 
30 
40 
50 
60 
7C 
80 
90 
10G 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
27C 
280 
290 
300 
310 
320 
33C 
340 
350 
360 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
52  0 
530 
540 
550 
560 
570 
580 
590 
600 
610 
620 
630 
640 
650 
66  0 
670 
680 
690 
700 
710 
720 
73C 
740 
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1000 


1100 


1200 


13  00 


1400 
1500 
1600 


1700 


$ JiSSiS 

$ 

$      THIS 

$AL0NG 

S 

kount  = 

NUMBER 
EPS   =    1 
TJ-ETA    = 
Y   =    YNE 
J    =    1 
XHDLO    = 
YhOLD    = 
UHOLO    = 
VhOLO    = 
MPMMM    = 
J   =    J+l 
CALL    XM 
CALL    VA 
IF    (J.E 
If    (J  .G 
IF    <EXA 
V0L3    = 
ThETAO 
CALL    RO 
GO   TO     1 
ThCNE    = 
THTWO    = 
VONE    = 
VTWO    = 
THNEW    = 
CALL    XM 
CALL    VA 
IF    (OAB 
CHECK    = 
IF    (CHE 
THTWO    = 
VTWO    = 
GG   TO    1 
ThONE    = 
VCNE    = 
CCNTI  NU 
GO   TO     1 
ACO    =     ( 
AOOA    = 
SUM    =     S 
SUMA    = 
IF    (KKK 


i$i$$$s$i$$$$s$iisi$$ti$t$iist$$$ss$ss$is%$$$$$$s$s$$ss$i$ 


SEGMENT    OF    THE     PROGRAM    PERFORMS    MjM 
THE    V    =    0    CONTOUR. 


ER1CAL    INTEGRATICN 

LEGtND 


=    5 

.0-11— 

c.o+o 


a 


X 

Y 

UNEW 
VNEW 
0 

ARCF    (XHOlD.YHOLD, THETA,  X,  YJ 

LUE     (X,Y,T  ,AL,U,V,MMMMMJ 

Q.2)     SIGN    =    V/DASS(V) 

T  .2  J    EXAM    =    V*SIGN 

M.LE.O.Q+G)    GQ    TC    1200 

V 

=    THETA 

TATE    (LeV£L,SIGN,TheTA) 

ICO 

THETAO 

THETA 
VOLO 
V 

({ (THONE-THTWCJ/JVTWQ-VONEJ )*VONE )+THCNE 
ARCH     iXriJLOf  YHOLOf  THNEWf  XfY) 
LUE     (X,Y ,T ,AL,UNEWf VNEW ,MMMMM) 
S(  VNEWJ.LT.l.O-liJ     GQ    TO    1600 

VNEW*SIGN 
CK.GT. 0.0+01    GO    TO    1400 

THNEW 
VNEW 
500 

THNEW 
VNEW 
E 

200 

(UNEW+UHJLOi*(Y-YhOLO)+(VNEW+VHOLO)*(X-XHOLD)  J  /TWOPI 
OABS( AUG) 
UM+ADD 
SUMA+AOOA 
KK.EQ.l)     WRITE     (6,1800)     ADC, AOOA 


**44******** *********** 

4 
IT    OF     THE    ALGCRI  THM.  * 


%\*44**4  ********  4**4*44*  *****************  Hit*' 
• 

*  THIS     IS    THE    TERMINATION    CRITERION    SEGMEN' 
I 

TERM    =    OABS(ADO) 
IF    (TERM.LT.E°S)    K0UNT  =  K.0UNT*1 

IF    (TERM.oc.EPSJ     KQUNT=0 

IF     (KOUNT.  EQ.NUMBEh)      GO    TO    1700 

*  « 

4*4************444*44*4444**44****4**4*****4*4**4****4*4**44****** 

ThETA    =    THNEW 

EXAM    =    -EXAM 

60    TO     1000 

WRITE     (6,2000) 

WRITE     (6,2100) 

CALL    GETi^E     ( I  ET  ) 

EL    =    DFLGAT(IcT)*2.6C-5 

WRITE     (0,1900)     EL 

*  $ 

STOP 


SUM, SUMA 

MMMMM 


I 

75  0 

76  0 
770 
780 
790 
800 
810 
820 
830 
840 
850 
860 
870 
880 
890 
900 
910 
920 
930 
940 
950 
960 
970 
98  0 
990 

ICOO 
1010 
102  0 
1030 
1C40 
1050 
1C60 
1070 
1080 
If  70 
1100 
1110 
1120 
1130 
1140 
1150 
1160 
1170 
1180 
1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1290 
1300 
1310 
1320 
1330 
1340 
1350 
1360 
1370 
1380 
1390 
1400 
1410 
1420 
1430 
1440 
1450 
1460 
1470 
1480 
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18C0 
1900 
2000 
2100 


FORMAT    (  1X.2E30. 10,/J 

FORMAT     (IX,  'ELAPSED    TIME 

FCRMAT     (lX,2E20.ol 

FORMAT    UXtiaj 

END 

SUBROUTINE     XMARCH     (XHCLC  ,YHOLD ,THET A,X , Y 1 

IMPLICIT    REAL*8    (A-H.O-Z) 

X    =    XH0L0-2.5l)-3*CC0S  IThETA) 

\T  =    YHOLD+2.  50-3*0  SI  NITHETA1 

RETURN 

ENO 

SUBROUTINE     ROTATE    (LEVEL ,SIGNf ThETA J 

IMPLICIT    °EAL*8    lA-H,0-Z) 

IF    (LEVEL. NE-1J    oO   TO     100 


DURING    THE    COMPUTATION    IS',E20.5) 

LEG  EN 


5 
5 


100 
2  00 

2. 

4. 
5. 


IF    (  SIGN. LT. 0.0  +  0) 
IF    (SIGN.GT.O.O+OJ 

GC   TO   200 

IF    ( SIGN. LT. 0.0+01 

IF    (S  IGN.GT  .0.0+0  ) 

RETURN 

END 


THETA=THETA-5.C-2 
THETA=TH£TA+5.D-2 

THETA=THETA+5. C-2 
THETA  =  ThETA-5.C-2 


LEGEND 


1490 

1500 

1510 

1520 

1530 

1540 

10 

20 

30 

40 

50 

60 

10 

20 

30 

40 

50 

60 

70 

80 

90 

100 


Mandatory  user  output. 

May  be  altered  by  the  user  at  his  discretion. 

Polar  radius  R  =  .025  built  into  program,  may 
be  changed  by  user. 


80 


SECTION  A-3.   SHANKS'  ACCELERATOR 
USER  INSTRUCTIONS 


The  algorithm  of  this  section  is  very  similar  to  that 
of  Section  A-l.   Only  those  features  which  differ  have  been 
discussed  in  this  instruction. 

1.  Contour  integration  is  performed  using  a  simple  para- 
meterized curve  for  a  finite  number  of  points  on  this  con- 
tour.  This  finite  number  of  terms  must  be  prescribed  by 
the  user.   The  integer  variable  M  is  used  for  this  purpose. 

2.  The  dimension  statements  at  the  beginning  of  the  program 
must  provide  storage  space  in  linear  arrays  which  are  equal 
to  or  larger  than  the  value  assigned  to  M.   (These  quanti- 
ties will  subsequently  be  used  with  Shanks'  accelerator. 
Note  particularly  that  the  array  is  doubly  indexed.) 

3.  The  integer  variable  N  specifies  the  number  of  iterates 
of  Shanks'  accelerator  which  are  desired  by  the  user.   How- 
ever, the  assigned  value  of  N  cannot  be  greater  than  6  as 
the  program  is  now  written. 

4.  A  sample  program  output  appears  in  Chapter  VI  of  this 
thesis.   Column  1  represents  the  natural  result  obtained 
after  M  evaluations  of  the  addend,  without  application  of 
the  accelerator. 
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I 

C              THIS    PROGRAM   CALCULATES    THE    INVERSE    LAPLACE    TRANSFORM    OF    A  10 

C              FUNCTION   USING    A    FINITE    NUM3ER    OF    TER'4S    ALONG    A    SIMPLE  20 

C               PARAMETERIZED    CONTOUR    CF    INTEGRATION    WITH    SHANKS'    ACCELERATOR  30 

IMPLICIT    '^EAL*d    (A-rt,0-Z)  40 

DIMENSION    AllOtlOOJ     f     ACD(IOO)  50 

DIMENSION    XXC100)  tYYClOOi  tUUJIOOJ  fVVdOOJ  60 

CALL    ERRSET    (207,256,-1,1,1.209)  70 

MCMMM    =0  .80 

TfcOPI     =■    4.0+0*0ARSIN(  1.0+0)  90 

AL   =    1.250-1  100 

T   =    TWOPI  110 

AA   =    3.0-1  120 

P    =   0.0+0  130 

M    =   20  140 

NN   =    0  150 

N    =    5  160 

CALL    CURVE     (AA, P,XO, YOJ  170 

CALL    VALUE     ( XO, YJ  ,  T , AL , LO , VO , MMMMN }  180 

C  190 

C                 $$$$SSSiS$$SS$i$iiS$$<i<<J$$SS$«S<$i$i5$SJS$S$$$$$l$<>i$S$JS$$S5v$5  20  0 

C              $  $              210 

C  220 

C               S      THE    FOLLOWING    SEGMENT    OF   THE    PROGRAM    PERFORMS    THE    CONTOUR  $               230 

C               SINTEGRATION    FOR    A    FINITE    NUMBER    CF    TERMS    USING    THE    SIMPLE  $               240 

C              iPARAMETERIZEO    CURVE.  $               250 

C               $  S               260 

200    NN   =    NN  +  1  270 

CALL  INCPEP  <P)  280 

CALL  CURVE  (AA,P,X,Y)  290 

CALL  VALUE  (A,Y,T,AL,U,V,MMMMMJ  300 

ACD(NN)  =  (  (V+VO)  *(X-XO)*(U+JO  )*(  Y- YO ) ) /TWOPI  310 

XO  =  X  320 

YO  =  Y  330 

UC  =  U  340 

VO  =  V  35  • 

XX(NN)  =  X  360 

YY(NN)  =  Y  370 

UU(NN)  =  U  380 

VV(NN)  =  V  390 

IF  (NN.EQ.M)  GO  TO  300  400 

GO  TO  200  410 

300  CONTINUE  420 

C               i  $      430 

C                 $$$$$$$S$Si$Siii$i$SS$$SS$iSSS$$Sli5iS*iii<i$iiiiii<J$$I$i$S$Jiii,£  440 

C  450 

Q                     *****«**«*4**#»**4*»**»»»***«*******«***********»*»*********»*****  A60 

C               *  *                470 

C               *      THIS    SEGMENT    OF    THE    PROGRAM     IMPLEMENTS    SHANKS'     ACCELERATOR    IN    *  480 

C              *ACCORDANCE    WITH    EQUATION    28    OF    CHAPTER    fc   OF    THIS    THESIS.  *               490 

C               *  *               500 

8   =    1.0+0  510 

EFS    =    1  .0-12  520 

C  530 

C  540 

DC   400    1=1, N  550 

C  560 

C  570 

DC  400    J=l,M  580 

400   A  (I,  J)     =    0.0+0  590 

C  600 

C  610 

1*1  620 

C  63  0 

C  640 

t50 
660 
670 
680 
690 

500   CCNTIMUE  700 

C  710 

C  720 

C  73C 

C  740 


DO    500    J=l ,M 

JM   =•  J-l 

8    =    ADO(J) 

IF    (  J.EC.l)     Ad  ,1 J 

=  e 

IF    (J.GT.l)    A< 1, J  ) 

=    8+A(l,JMJ 

CCNTIMUE 
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00  600    1=2, N  750 

IM   =    1-1  760 

K    =    2*1-1  770 

C  780 

C  790 

DO   600    J=K,M  800 

Jf   =    J-l  810 

JMM   =    J-2  820 

A<1, J)     =    l.D+O  e30 

DEN   =    A(IM,  Ji  +  A(I  M,  JMM)-2.C+0*A(  If,  JMJ  840 

IF    (OABS(DEN) .GT.EPS)     A ( I , J) =( A( I  V  ,  J) *A ( I C , JMf ) -A ( 1 f, JM )**2) /OEN  850 

600    CCNTINUE  860 

C  870 

C  880 

C  890 

C  900 

OC    700    J=2,M                                                                              .  910 

700    WRITE     (6,800)    J , < Al I , J  J , 1=1 ,N )  920 

C  93  0 

C              *                                                                                                                                                                            *  940 

C  960 

STOP  970 

C  980 

C  990 

800    FCRMAT     <5X  ,  13,6 Fl 8  .14  J  1C00 

ENO  1010 

SUBROUTINE    CURVE    (A,P,X,VJ  10 

IMPLICIT    REAl*8     (A-H.0-2J  20 

X   =    A-(P*P)  30 

Y    =    P  40 

RETURN  50 

END  60 

SUBROUTINE     INCREP     (P)  10 

IMPLICIT    REAL*8    IA-H.C-Z)  20 

P    =    P  +  l.D-1  30 

RETURN  40 

END  50 
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SECTION  A-4.   SPECIAL  CONTOUR  ALGORITHM  FOR  CASES  AS 
TREATED  IN  CHAPTER  EIGHT 

USER  INSTRUCTIONS 


The  integration  contour  used  in  this  section  is  a  vertical 
path  to  a  specified  value  of  y,  followed  by  a  simple  parabola. 
The  algorithm  of  this  section  is  very  similar  to  that  of 
Section  A-l.   Only  those  items  which  differ  are  discussed  in 
this  section. 

1.  Contour  integration  is  performed  in  the  initial  segment 
of  the  program  along  the  Bromwich  contour.   The  integer  vari- 
able JJJJJ  specifies  the  number  of  poles  of  the  function 

ct  s 
f(s)  =  —  tanh(-?r)  above  the  real  axis  that  are  to  be  enclosed 
v      s       2 

to  the  left  of  the  distorted  Bromwich  contour,  which  is  used 
in  the  second  segment  of  the  program. 

2.  The  variable  HALT  computes  the  elevation  above  the  real 
axis  corresponding  to  the  input  value  of  JJJJJ  for  the 
specified  hyperbolic  tangent. 

3.  The  variable  FLAG  is  internally  assigned  a  value  of  zero 
or  one,  corresponding  to  the  first  or  second  segments  of  the 
contour  of  integration  which  the  algorithm  employs.   No 
action  is  required  by  the  user. 

4.  The  integer  variable  KSTART  is  utilized  in  the  outer  loop 
of  the  program  to  increment  the  value  of  t  for  which  F(t) 

is  computed. 

5.  The  linear  arrays  STORE  and  STASH  are  used  to  convert  the 
values  of   F(t)  and  t  to  single  precision,  which  is  required 
for  the  graphics  package  associated  with  the  IBM  360/67. 
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6.   SUBROUTINE  VALUE  contains  the  real  arithmetic  coding  of 
the  function  g(s)  for  the  hyperbolic  tangent  function  of 
Chapter  VIII  of  this  thesis.   Standard  trigonometric  identi- 
ties for  the  hyperbolic  functions  of  a  complex  argument  have 
been  utilized. 

Note:   This  procedure  is  obviously  keyed  to  the  particular 
function  f(s)  =  —  tanh(^).   The  user  must  modify  it  for 
other  functions.   One  obvious  modification  which  would 
generalize  the  procedure  would  simply  be  to  specify  rather 
than  calculate  the  value  of  HALT,  i.e.,  input  HALT  rather 
than  JJJJJ. 
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C  THIS    PROGRAM   COMPUTES    THE    INVERSE    LAPLACE    TRANSFORM    OF    A 

FUNCTION    USINu    THE    SIMPLE    PARAME TERI ZEO    CONTOUR    OF     INTEGRATION 

C  WITH    MODIFICATION    TO    VAPY    THE    tLcVATION    FROM    WHICH    THE    CURVE     IS 

C  LAUNCHED.        THE    BRGMWICH    CONTOUR     IS    USED     IN    THE     INITIAL    SEGMENT. 

INPLICIT    RE;L*8     (A-h, G-Z)  20 

REAL*4    STORE,     STASH  30 
DIMENSION    STORE! 50) t     STASH! 50) 

CALL    ERRSET     (z07 , 256 , -i  ,  1, 1 ,209  )  50 

TWOPI     =    4.D+0*0ARSINI  1.0  +  0)  60 

AL   =     1  .0  +  1  70 

JJJJJ    =    15  80 

KKKiU    =0  90 
NUMBER    =    5 
KSTART     =   50 

C  100 
00    400    KKK    =     1,     KSTART 

FLAG    =    0.0+0  120 

HALT    =     (OFLOAT(  JJJJJ)  *Ti»0PI)/(2.D+0*AL)  130 

ACO   =    0  .0+0                                                                                                                                   -  140 

ACOA    =    0.0+0  150 

SLM    =0.0+0  160 

SUMA    =    0.0+0  170 

LCSC    =    0  180 

MMMMM    =    0  190 

AK    *    KKK  20C 

AA    =    3.0-1  210 

T    =    5.0-1+AK*1.0+0  220 

P    =    0.0+0  230 

KCUNT    =    G  240 

NN   =    0  250 

EPS    =    1  .D-ll  260 

OELP   =     1.0-1  280 

CALL    CURVE     ( AA, P , HALT , F LAG , XO ,Y0)  290 

CALL    VALUE     (XO  »Y0  tT  tALtUC»VOf  MMMMM)  200 

TIME    =     C.0+0  320 

CALL    SETIME  330 

2  00    NN    =    NN+1  340 

CALL    INCREP     (  HALT  ,06  LP  ,  F  .FLAG  J  350 

CALL    CUPV6     (AA,P,  HALT  ,  FlAG,X,  Y)  360 

CALL    VALUE     IXfYfTfALtU»ViMMMMMJ  270 

ADO   =     (  (V+VO)  *(  X-XJ  )  +  (  U  +  UO)  *(  Y-YO  ))  /TWOPI  380 

IF    (NN.EC.l)    AOJOLO=A0C  39C 

PROO    =     AOG*ADOOLJ  400 

IF    (PROO.LT.G.J+0)    L0SC  =  L0SC+1  410 

AOOA    =     OABS (AOO)  420 

SOM    =     SLM+AOD  430 

SUMA    =     SUMA+AOOA  440 

TEST    =    OABSIAJO)  450 

IF    (TEST.LT.EPS)     KOUNT  =  KCUNT+ 1  460 

IF    (TEST.GE.EPS)     KJJNT=0  470 

CALL    GETIME    (iET)  490 

EL   =    DFLOATd  ET)*2.6D-5  500 

TIME    =    TIME+EL  510 

CALL    SETIME  520 
IF(KOUNT.£Q. NUMBER)     GO    TO    300 

XC   =    X  540 

YO   =     Y  550 

UO   =    U  560 

VC    =    V  57G 

ADOOLO    =    AOO  580 

GO    TO     200  590 

300    CCNT  INUE  600 

WRITE     (6,500)    NN,TI*E  610 
WRITE     (6,600)     T,     SUM,     SUMA,     X,     Y,     LOSC,     MMMMM 

STORE(KKK)     a    SNGL(T)  630 

STASH  (KKK)     =    SNGL(SUM)  640 

400    CONTINUE  650 
CALL     PL  CTG  (ST  ORE,  ST  AS  (-.50, 1,1,0,  • T«,i,«F(T)  ',4,0.0,  55. 0,-2.0, 2. C, 
16.0,  6.  0) 
CALL    PLOT  (0.0,0.0,999) 

C  660 

STOP  670 

C  680 

500  FORMAT  (IX,  'TOTAL  ELAPSEO  TIME  AFTER  M5,'  ITERATIONS  IS  ',£20. 5)  690 
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6  00 


100 

2  00 


FORMAT 
ENO 

SIBR3UT 
IPPLICI 
IF  (P.L 
X    =    AA- 

Y  =    P 
GC   TO    2 
X    =    AA 

Y  =     P 
RETURN 
ENO 

SUBROUT 
IPPLICI 
P  =  P  +  0 
IF  (P.E 
RETURN 
ENO 

SUBROUT 
IKPLICI 
HMMMM    = 
TEM1    = 
TEM2    = 
TEM3    = 
TEM4    = 
TEM5    = 
TEM6    = 
TEM7    = 
TEM8    = 
TEM9    = 
TEM10    = 
TEM11    = 
TEM12    = 
TEM13    = 
TEM14    = 
U    =    (TE 

Y  =    <  TE 
RETURN 
ENO 


(IX,  1P5E20.5.2I8,/) 

INE    CURVE    (AA, F ,HALT,F L AG , X  ,Y) 

T    REAL*3     (A-h,0-Z) 

T. HALT. AND. FLAG. NE. 1.0+0)     GO    TO    100 

( (P-hALT )*(P-HALT)  ) 

00 


INE     INCREP     (  HALT, OELP,P, FLAG) 

T    REAL*3     (A-H,C-2J 

ELP 

Q.HALT)    FLAG=1.0+0 


INE    VALUE     (X, Y,T,AL,U,  V.MMMMM) 
T    REAL*3     IA-HtG-ZJ 

MMMMM-t-1 
(AL*X) 
(AL*Y) 
DSINH(TEMl) 
DSIN(TEM2) 

0CCSH(TdMl)*DCCS(TEM2) 
TEM3/TEH5 
TEMWTEMt? 
(TE^6*x)+(  TEM7*Y) 
(TEM7*X)-<  TEM6*Y) 

(X*X)  +  (  i*1i 

TEM8/TE.1ia 

TEM9/TEM10 

PCXP(  X*TJ*OCCS( Y*T) 

LCXP(  X*T)  *OSIN(  Y*T) 
M13*TEMU1-<TEML4*TEM12J 
M14*TEM11)  «■(  TEM13*TEM12) 
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APPENDIX  B 
LISTING  OF  TEST  CASES  AND  EXPERIMENTAL  RESULTS 

The  Laplace  transformation  pairs  f(s)  and  F(t)  which  we 
have  used  as  test  cases  are  listed  within  this  appendix. 
Also  included  for  each  function  f(s)  is  the  listing  of  the 
SUBROUTINE  VALUE  and  a  typical  result  of  numerical  inversion 
at  one  specific  value  of  t.   Case  9  also  shows  a  plot  of 
F(t)  versus  t  obtained  numerically. 

The  results  obtained  by  numerical  inversion  have  been 
compared  where  applicable,  with  tabulated  values  listed  in 
(1).   In  all  other  cases  the  value  of  F(t)  has  been  cal- 
culated on  the  IBM  360/67. 

The  symbols  which  are  used  in  the  presentation  of 
results  are  defined  below: 
AL    -  constant  associated  with  f(s) 
BE    -  constant  associated  with  f(s) 
T     -  time  for  which  F(t)  is  evaluated 

AA    -  real  axis  starting  position  of  the  integration  contour 
x     -  final  value  of  x  at  termination 
yf    -  final  value  of  y  at  termination 

MMMMM  -  total  number  of  calls  made  to  SUBROUTINE  VALUE 
LOSC   -  number  of  sign  changes  between  successive  addends 
SUM   -  value  of  F(t)   (equation  (16)) 
SUMA  -  value  of  G(t)   (equation  (19)) 

TIME   -  clock  time  (sec)  required  for  numerical  integration 
ANALYTICAL  -  analytical  value  of  F(t) 
In  all  cases  we  used  N  =  5  and  EPS  =  l.D-11. 
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Case    1 . 


X 


-1 


SU3R0UT 
THIS     IS 
1MPLICI 
I  MPL I  C  I 
MMMMM    = 
XPTR    = 
XPTI    * 
CXPT    = 
CNUM   ■ 
CDEM3M 
CRESLT 
U   =    ORE 
V    =    DIM 
RETJRN 
ENO 


INE    VA 

TABLE 

-     Q  P5L 

T  C3MP 

X*T 
Y*T 
DCM'LX 

CQEXPI 
=  DC  MP 
s  C^'JM 
AL(C<?E 
AGCRE 


L'JE    (X,Y,T,ALtU,V,MMmmm) 

EiTRY    NUM3E-,     01    OF     LAPLACE 
*8     (  A,  3,  D-H,0-Z) 
LEX*lo    (C) 

+    1 


(XPT^,XPTJ    ) 
CXPT) 

LX( a,Y) 

/CDE'IOM 

SLT> 

SLT) 


TRANSFORMS    BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.2 

XF   =   -4.64 

y?  =  2.2 

MMMMM   =23 

LOSC    =    5 

SUM   =    1.00000D+0 

SUMA   =    1.42588D+0 

TIME   =   0.79794D-1 

ANALYTICAL   =    1.00000D+0 


93 


Case   2. 


X 


-I 


SUBOiUT 
THIS     IS 

i  mp  li  :  i 

IMPLICI 

MMMMM    = 
CS1    =    0 
XPTR.    = 
XPTI    = 
CX^T    = 
CNUM    = 
CDENOM 
CDEN^M 
CRESLT 
U   =    &AE 
V    =    DIM 
RETURN 
END 


IME 
TAB 

T    RE 
T    CO 

CMDL 
X*T 

DCM=> 

COEX 

=  oc 
=  cs 

=     CN 
AL(C 


VALUE     IXf YtT.ALtUtVfMMKMMI 

LE    ENTRY     N'JM3£R    02    OF     LAotACE    TRANSFORMS 
AL*8    ( AtBtD-HtO-Z) 
MPLEX*16    (  C  J 

MM     +     I 

x(  x,Y) 


8Y    SPIEGEL 


LX(X?TR,XPTI  ) 
3(CX^T  ) 
MPLX (X  ,Y) 
1*CS1 
'JM/CDtNOM 
RESLTJ 
AG  C  CRESLT) 


AL  =    1.25D-1 
T   =   8.D+0 
AA    =    0.2 
Xp   =   -4.6U 

yF  =  2.2 

MMMMM   =   23 

LOSC    =   k 

SUM   =   8.00000D+0 

SUMA   =   8.04195D+0 

TIME   =    0.532^8D-1 

ANALYTICAL   =   8.00000D+0 
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Case   3 


£ 


-1 


■±-  n  =  1,2.3, 


—  (n-D! 


0!  =  1 


SUBROU 

THIS  I 
IMDLic 
IVPLIC 
MM.MMM 
ZERO  = 
CS1  = 
XPTR  = 
XPTI  = 
CXPT  = 
CMUM  = 
CDENOM 
CRESLT 

u  =  on 

V  =  DI 
RETURM 
ENO 


TIME    VALUE(X,Y,T,il,NDlWEStU»V,M'4MVM) 

S    T^<?LE    ENTRY    NUMBER    03    OF     LAPLACE    TRANSFORMS     8Y    SPIEGEL 

IT   o  <EiL*6    (  A,3.D-H,0-Z) 

IT    COMPL£X*i6     ( C) 

=    MMMMM    ♦    1 

0.0+0 
0CMOLX(X,Y) 

X#T 

y*r 

DCMPLX  (XPTR.XPTI) 

COEXP(CX»T) 

=    CSl**N°nwER 

=    CNUM/CDtNJM 
E4L(CRE5LT) 
MAG(CRE5LT) 


AL   =    1.25D-1 
N POWER   =    3 
T   =   8.D+0 
AA   =    0.3 

X      =    -4.54       , 

r 

JF 

MMMMM   =   23 

LOSC    =    3 

SUM   =    3.20000D+1 

SUMA   =    3.20008D+1 

TIME  =   0.8369^0-1 

ANALYTICAL   =    3.2000CD+1 
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Case  k. 


£ 


-1 


n  >  0 


T(n) 


SUBROUTINE    VALUE     (X,    Y,    T,    &L, 

POWED,    Ut 
OF    LAPLACE 

THIS    IS    "M8LE    ENTRY    NUMBER    04 

IVPLICIT     REAL*8     (  A,  3,  C-t-,l-Z  ) 

I  MPLICIT    C0MPLEX*i6     <C} 

DIMENSION    MCALLC5) 

MMMMM      =      MMMMM     +     1 

CS1    =    DCMPLX(X,YJ 

XPT»    =     X*T 

XFTI     =    Y*T 

CXPT    =    DCMPLX(XPTR,XPTI) 

CMUM    =     CDEXP(CX?T ) 

CALL   CPCWER    (CSltCOENCM, POWER, 

l.NCALL) 

CRESLT    =   CNUM/C3EN0M 

U    ■    DREAL(CPESLT) 

V   =    DIMAG(CRESLTj 

RETURN 

END 

V,    j|mmmm,    NCALL) 
TRANSFORMS    eY    SPIEGEL 


AL  -    1.2  5D-1 
POWER   =    1.7D+0 
T  =   8.D+0 
AA   =    0.2 
XF   =   -0.6U 

y«  =  2.2 

MMMMM   =23 

LOSC    =   I* 

SUM  =  U.71815D+0 

SUMA  =  U.832UD+0 

TIME  =  0.85098D-1 

ANALYTICAL  =  J+.71815D+0 
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Case   5 


JS 


-1 


««< 


SUBROUT 
THIS  IS 
IMPLICI 

impli:i 

HMMMM     = 

ZERO    = 
CSl    =    D 
XPTR    = 
XPTI     = 
CXPT    = 
CMJM   = 
CS2 
CDENOM 
CRESLT 
U   =    ORE 
V    =    DIM 
RETJRN 
END 


=    D 


IME     VALUE     (X.Y.TjftL.'J.V.MMWVMJ 

T«?LE    E^JTRY    NUMBER    05    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 
T    REAL*8     IA,BtD-H,0-ZJ 

T   COMPLEX*  16     (C) 

MMMMM      +     I 

0.0+0 

CM»LX(X,Y) 

X*T 

y*t 

DCMPLXCXPTR.XPTI) 

CDEXP(CXPT  ) 
CMPLX(ALt ZERO) 
=    CSl   -    CS2 
=    C'lUM/COENOM 
AL(CRESLT) 
fiG(CPESLT) 


AL   =    1.25D-1 
T   ■   8.D+0 
AA    =    0.3 
XF   =   .^. 99 

yF  =   2.3     2  '"* 

MMMMM   =   23 

LOSC    ■   5 

SUM  =  2.71828D+0 

SUMA   =   3.59U74D+0 

TIME  =   0.89050D+0 

ANALYTICAL   =   2.71828D+0 
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Case  6 


X 


-1 


(*-«)» 


n  =  1.2,3, 


(n-1)!' 


0!  =  1 


SUBROUTINE    VALUE(X,Y,-, 4L»NPOWER,UtV, MMMMM J 

THIS    IS    TABLE    E^ITRY    NUMBER    06    IF     L4°LACE    TRANSFORMS 

IMPLICIT    REAL*8     (  A,  Bt  D-l-,0-Z  ) 

IMPLICIT    C0"PLEX*16     (CJ 

M^IMMM    =     MMVJMM     +     i 

ZERO    =    0.0  +  0 

CS1    =    0CMOLX(  X,Y) 

CS2    =     OCMPLX(AL,ZERO) 

XPTP.    =    X*' 

XPTI     =     Y*T 

CXPT    =     DCMPLX(XPTR,XPTI ) 

CNUM   =    CDEXP(CXPT) 

CS3   =    CS1   -    CS2 

CDENOM    =    CS3**NP0WER 

CRESLT    =   CNUM/COENOM 

U    =   OREAL(CRESLT) 

V    =    OIMAG(CRESLT) 

RETURN 

END 


BY    SPIEGEL 


AL   =    1.25D-1 
NPOWER    =    3 
T   =    8.D+0 
AA   =    O.U 
Xp  = -4.44 
yF  -  2.2 

MMMMM   =   23 

LOSC    =   3 

SUM  =  8.69850D+1 

SUMA  =  8.69365D+1 

TIME  =  0.75790D-1 

ANALYTICAL  =  8.6985OD+I 


98 


Case   7 


£ 


-1 


(s-a)" 


n  >  0 


f-1  e"' 


SU8R01' 
THIS  I 
IMPLIC 
I ^PLIC 
DIMEMS 

Mr""" 

ZERO   = 
CS1    = 
CS2    = 
CS3    = 

XPT5 

XPTI  = 
CXPT  = 
CNM  = 
CALL  C 
CPESLT 
ij  =  OR 
V  =  01 
RETJ'M 
END 


TINE 
S    T4 

IT  R 
IT  C 
ION 

=      V|M 

C«  D 

DCMD 

DCMP 

:si 

X* 

Y* 

DCM 

COE 
PCWE 

=  r 
EAL( 
"AG< 


VAL'J 
BLE  E 
E£L*8 
HVPLE 
VC  4LL 
"MM  + 
+  0 

LX(  X, 
LX  (AL 
-    CS2 


c     (X ,Y  t TtALrNPOWER ,U ,V , MMM"Mt NC4LL ) 
NTRY    NUMBER    07    CF    LAPLACE    TRANSFORMS 

(  A,?,  C-h,0-Z) 
X*16     (C) 
(5) 

1 

Y) 

,ZER0J 


BY    SPIEGEL 


PLX( XPTRi XPTI ) 

X  P  (  CX  °  T  J 

R     (CS3,  COENCM,  POWER,  1  ,NCALU 

MiJM/CDEWJM 

CRESLTJ 

CRESLTJ 


AL   =    1.25D-1 
POWER   =    1.7D+0 
T   =   8.D+0 
AA   =   0.3 

x  =-4.54-  ~t-^u 

F 
yF   -  2.2 

MMMMM  =  23 

LOSC  =  k 

SUM  =  1.28253D+1 

SUMA  =  1.30533D+1 

TIME  =  0.96928D-1 

ANALYTICAL  =  1.28253D+1 
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Case   8. 


£ 


-1 


sin  at 


■»-»■  a* 


SU8R1 
THIS 
IMPLI 
IfOLI 

MMMMM 

ZER3» 

CS=DC 

CALO 

CT=0C 

CDEN 

CEXP  = 

C    =    C 

U=OPE 

V=DI^ 

RETJR 

END 


UT  IN  E 
I  S  XA 
CIT  * 
CIT    C 

=      MM 

0.0*0 
MPLX( 
CMPLX 
^PLX( 
=     (<C 

:dex» 

EX  P/C 
ALIO 
fiG(C) 
N 


VALJEt  X,Y, T, 41, U, V,MMMMM) 
BLE    E.HTW     'vlUMEcS    03    OF     LAPLACE    TRANSFORMS     BY    SPIEGEL 
EAL*8    ( A, 6  ,0-H,0-Z) 
?'1DLEX*I6     (C) 

MMM     +     I 

X,Y) 

(  *L,ZE*0) 

T,ZER0) 

C*CS)-M  CAL*CAD) 

(CS*CT) 

DEN 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.3 
xF   =   -4.54 
yF  -  2.2 

MMMMM   =   23 

L05C    =   4 

SUM   =   6.73177D+0 

SUMA  =  7.03097D+0 

TIME  =  0.86294D-1 

ANALYTICAL  =  6.73177D+0 
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Case   9. 


£ 


"1 


*2  +  a2 


cos  at 


SUBR^UT 
THIS     IS 

impli:  I 

I^IPLICI 

MMMMM  = 
ZEPD  =  0. 
CS=OCMP 
CAL=DCM 

cto:mp 

CDEN=C3 
C=CS/CO 
CExt>=ro 

C=C*CEX 

U=D*E'-L 

V=OHAG 

RETURN 

END 


IME 

T   5 

T 

D+0 

LX( 

PLX 

LX( 

**2 

EN 

EXP 

P 

(C) 

(C) 


V»LJE     ( Xf     Y,     T,    1L.     U,    V,    MMMMM) 
BLE    ENTRY    MUM  BE0    09    ?F    LAPLACE    TRANSFORMS    BY    SPIEGEL 
HAL*8     ( A,o,0-H,a-Z) 
IMPLEX* lo    ( C  ) 
MMM    +    1  * 

X,Y) 

{ AL,ZERG) 

T,ZERD) 

+  CAL**2  , 

(CS*CT) 


AL   =    1.25D-1 
T   =    8.D+0 
AA   ■    0.3 

xr  =  -*.99 

yF  -  2.3 

MMMMM   =   2b 
Lose  =  5 

SUM   ■    5.^0302D-1 
SUMA   ■   2.00010D+0 
TIME   =    0.11138D+0 
ANALYTICAL   =    5-^0302D-l 


30. 

T 
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Case    10. 


£ 


'1 


1 

(»  -  b)-  +  a- 


e'-"  sin  at 


S'JBROU 

THIS    I 
I"-1PLIC 
I VPLIC 
MMMMM 
ZERO    = 
CS1    = 
CS2    = 
CS3    = 
XPTR    = 
XPTI     = 
CXPT    = 
CNUM    = 
CS4    = 
CDENOM 
CR6SLT 
U    =    OP 
V    =    01 
RETURN 
END 


T IME    VALUE(X,YtT,ALt3E,UfVf MMMMM) 

S    TABLE    ENTRY    NUMBER    10    OF     LAPLACE    TRANSFORMS     8Y    SPIEGEL 
IT    R  EtL-8    (  A,3,0-H,0-Z) 
IT    Cn^PLEX-li     Id 
=    MMMMM    +     1 

0.0  +  0 
DCMDLX (X,  Y) 
DCMPLX(BE,ZEROJ 
DCMPLX( AL  ,Z£R0  J 

X*~ 

Y*T 

DCMPLX(XPT^tXPTI ) 

C0EXP(CXP-  ) 

: si  -  csz 

=    (CS4*CS4.)  +  (CS3*CS3) 

=    CN-JM/COtNOM 
EAL(CRESLT) 

MAG(C°ESLT) 


AL   =    1.25D-1 
BE   =    1.25D-1 
T   =   8.D+0 
AA   =    0.3 

xF  =  -4.5* 

yF  =  2.2 

MMMMM   =23 

LOSC    =   k 

SUM  =   1.82988D+1 

3UMA   =   1.88117D+1 

TIME   =    0.700^4D-1 
ANALYTICAL    =    1.32988D+1 
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Case    1 


.r1 


*-  b 


(«-*)»  +  0* 


c6'  cos  at 


SUBROUT  INE    VALUE (  X,  Y,  T  ,AL,BE  ,  U,  V,  mmmm^) 

THIS     IS    TARLE    ENTRY    N>W«£R    11    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    REAL*8    ( A, 8 , J-H  ,0-Z) 

IMPLICIT    C-»^DLEX*lo    (  C  ) 

MPMMM    =     MMMMV    +    1 

ZERO    =     0.0+0 

CS1    =    DCM»LX(XfY) 

CS2    =    DCMPLX( ec,ZE"0) 

CS3    =    OCMPLXtALf  ZERO) 

XPTR    =    x*t 

XPTI     =    Y*T 

CXPT    =    DCMPLX(XPTR.XPTI) 

CMM    =     (CS1-CS2)*C0EXP(CXDTJ 

CS4   =    CS1    -    C52 

CDEN1M    =    (CS**CS4)+<CS3*CS31 

CRESL'     =    CMUM/CDENOM 

U  =    DREAL(CRESLT) 

V    =    DIMAG(CRESLTJ 

RETURN 

END 


AL   =    1.25D-1 
BE   =    1.25D-1 
T   =   8.D+0 
AA   =    O.U 
Xp    =    -4.89 

yF=    2.3 

MMMMM   =   24 

LOSC   =   5 

SUM  =    1.46869D+0 
SUMA   =  4.61766D+0 
TIME   =   0.78130D-1 
ANALYTICAL   =    1 . U6869D+O 
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Case    12. 


£ 


-1 


—  sinhnt 


SUBROUT 

THIS  IS 
IMPLICI 
I  VPUCI 
MMMMM  = 
2E"1=0. 
CS*OCMP 
CAL=l)CM 
C'=DCMP 
COEN    = 

r exp=:d 

C  =  CEX 
U=OREAL 
V=DIMAG 

RETJRN 
END 


IN6    VA 

TABLE 
T    =  EiL 

T  co«p 

D+O 

LX(X,Y 

PLX( AL 

LX(Tf Z 

UCS*C 

EXP(CS 

P/CDEN 

(c; 

(C) 


L'JE(X,Yf  T,  AL,  Ut  V.^mmmm  ) 

ENTRY    MUMeER    12     OP    LAPLACE    TRANSFORMS    BY    SPIEGEL 
*8    ( A,B,0-H,0-Z) 
LEX*  16     (CJ 

♦    1 

i 


ER1) 

S)-<  CAL*CAD) 

*CTJ 


AL   =    1.25D-1 
T   =   8.  D+O 
AA   =    0.3 

xF  =  -b.51* 

yF  =  2.2 

MMMMM   =    23 

LOSC    =   b 

SUM  =   9.14.0162D+0 

SUMA   =   9.6^RR9D+0 

TIME  =   0.10691D+0 

ANALYTICAL   =    9.U0162D+0 
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Case    13 


X 


-JL 


cosh  at 


SUBROU 

THIS  I 
IMPLIC 

I *PLIC 

MMMMM 

ZERO=0 

CS=DCM 

CAL=0C 

CT=OCV 

COEN    = 

c=cs/: 

CEXP=C 
C=C*C6 
U=DRE& 
V=DI^A 
RETJRN 
ENO 


TINE 

S  TA 
IT  5 
IT  C 
=  MM 
.0+0 
PLX( 
MPLX 
PLX( 
((C 
DEN 
OEX' 
XP 
L(C) 
G(C) 


VALJE<X,/,T,  AL,U,  V.mmmmm  ) 

6LE    EMTRY    NUMBER    13     OF    LAPLACE 
Eil*8    {  A,  3  ,0-ri,>Z) 
OMPLEX*lo     (CJ 
MMM    ♦    1 

X,Y) 

(  ftL.ZERO) 

T,ZERO) 

S*CSi-ICAL*CAL)J 

(  CS*CT) 


TRANSFORMS  BY  SPIEGEL 


AL  =  1.25D-1 
T  =  8.D+0 
AA  =  0.3 
Xp  =  -^.99 

yF  -  2.3 

MMMMM  =  2U 

LOSC  =  5 

SUM   =    I.5U308D+O 

SUMA   =   2.52i+01D+0 

TIME   =   0.85176D-1 

ANALYTICAL   =    I.5U3O8D+O 
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Case    14 


X 


-1 


(*  -  by*  -  at 


r"  sinhnf 


SUBROUT 

THIS  IS 
IMPLICI 
IfOLICI 

"•1MMMM    = 
ZERO    = 
CS1    =    D 
CS2    =    D 
CS3    =     D 
XPTR    = 
XPTI     = 
CXPT    = 
CMJM   = 
C  «4   =    C 
CDENOM 
CRESL" 
U   =    ORE 
V    =    DIM 
RETURN 
END 


INE    VALJE< X, Y,T, AL, BE, Ut  V,MMMMM) 

T48LE    ENTRY    NUMBER    14    OF    U  =  LACE    TRANSFORMS 
T    =  E4L*8    (  Af<3  ,D-rl,0-Z)  - 

T    CO^PLEX*lo     (C) 

O.D+O 

CMPLXU,  Y) 

CMPLXCBE  .ZERO) 

CM°LX< AL.ZERO) 

X*' 

Y*T 

CCMOLX(X°TR,XPTI  ) 

COEXD(CXPT) 

SI   -    CS2 

=    (CS4*CS4)-<CS3*CS3) 

=    CNUM /CDENOM 

AL  (C=ESLT) 

AG(CPESLT) 


BY    SPIEGEL 


AL   =    1.25D-1 
BE   =    1.25D-1 
T   =   8.D+0 
AA   =    0.5 
Xp    =   -4.79 

yF  =  2-3 

MMMMM   =   24 

LOSC    =   4 

SUM   =   2.55562D+1 

SUMA   =   2.81818D+1 

TIMS   =   0.7^9320-1 

ANALYTICAL   =   2.55562D+1 
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Case    15 


X 


-1 


-  b 


(s-b)- 


ebt  cosh  at 


SUBROUT INE    VALUE(X,Y,T,AL,*E,UiV, VMMMM) 

THIS    IS    TABLE    ENTRY    NUMBER    15    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    5EiL*3    [  Af3fD-H,Q-Z) 

IMPLICIT    COMPLEX*16     (CJ 

MMMMM     =     MMMM.M    ♦     1 

ZERO    =     O.O+O 

CS1    =    OCMPLX(X,Y) 

CS2    =    0CMPLX(3d  f  ZERO) 

CS3    =    OC^PLX(AL, ZERO) 

XPTR    =    X*T 

XPTI     =    Y*T 

CXPT    =     OC^PLX  (XPTR,  XPTI) 

CNUM    =     1CS1-CS2 J*COEXP(CXPT) 

rS4   =     C  SI    -    r  S  2 

CDENOM    =    (CS4*CS+)-(CS3*CS3) 

CRESLT    =   CNUM/CDENOM 

U    =    DP  EAL(C°ESLT) 

V    =    OIMAG(CRESLT) 

RETURN 

ENO 


AL   =    1.25D-1 
BE   =   1.25D-1 
T   -   8.D+0 
AA   =   0.5 
XF    =   -4.79 

yF  =  2.3 

MMMMM  =24 

LOSC    =   5 

SUM   =   4.19453D+0 

SUMA    =   9.63599D+0 

TIME  =  0.77740D-1 

ANALYTICAL   =   4. 19^53D+0 
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Case    16 


£ 


-1 


i 

(»-a)(»-  *) 


#  6 


—  e"'  -  f" 

6-a 


SUBROUTINE     VALUE  <  X,  Y,  T  ,  M  ,3E  ,U  ,  V  f  wvmmm) 

THIS    IS    Ti*LE    ENTRY    NUMBER    16    :>F     LAPLACE    TRANSFORMS    8Y    SPIEGEL 

IMPLICIT    PEAL*3     ( A, 5, C-h,3-ZJ 

IMPLICIT   COMPLEX*  16    (C) 

MfMMM    =     MMMMM     +     1 

ZERO    =    O.D+0 

CSl    =    OCMPLX(X.Y) 

CS2    =    DCMPLX (BE.ZERH) 

CS3   =    DCMPLXULtZERO) 


XPTR 

=     X*T 

XPTI 

=     Y*T 

CXPT 

=     DC^PLX(XPTR, 

<°T 

I  ) 

CNUM 

=     CDEXP(CX^TJ 

CS4    = 

■    CSl    - 

CS3 

C  S5    = 

■    CSl   - 

CS2 

CDEN?M    =    CS4*CS5 
CRESLT    =   CN'JM/COENOM 
U    =    DREALtCSE  SL  T) 
V    =    DIM AG(CRESLT) 
RETURN 
END 


AL   =    1.25D-1 
BE   =   -1.25D-1 
T   =    8.D+0 
AA   =   0.4 
Xy   =  -4.39 

yF  =  2.3 

MMMMM  =24 

LOSC    =   4 

SUM  =   9.40161D+0 
SUMA  =   1.07718D+1 
TIME   =   0.76648D-1 
ANALYTICAL   =    9.40161D+0 


108 


Case    17. 


X 


-1 


(s-a)(a-b) 


**  b 


—  bc<><  -  ac"' 

A  -  a 


SU8R0UT  INE    V ALJ 5( X,  Y  ,  ?  ,  AL , 3E, J  ,  V  ,«M'<M1) 

THIS    IS    TA9LE    ErtTRY    .NUMBER    17     OF     LAPLACE    TRANSFORMS 

IMPLICIT    ?EAL*3     {  A, 9  ,D-HfO-ZJ 

IMPLICIT    COMPLEX*lb     (C) 

MMMMM    =     MMMVM    *■      i 

ZERO    =    O.O+O 

CS1    =    DCMPLX<X,Y) 

CS2   =    DCM°LX(BE  ,ZERO) 

CS3    =    DC^PLX  (AL,  ZERT) 

XPTR    =    X*t 

XPTI     =    Y*T 

CXPT    =     OCMPLX(X°TR,XPTI) 

CNUM    =    CS1*CDEXP(CXPT) 

CS4   =    CS1   -    CS3 

CSS    =    CS1    -    CC2 

COENOM    =    C  S4*CS5 

CPESLT    =    C'UV/C DENIM 

U    =    OREAL(CPESL'p) 

V  =    DIMAGICRESLTJ 

RETJRN 

END 


BY    SPIEGEL 


AL   =    1.25D-1 
BE   =    -1.25D-1 
T   =   8.D+0 
AA   =    0.3 

X^,    =    -4.99 

r 

yF  =  2.3 

MMMMM   =24 
LOSC   =   5 

SUM   =    1.54308D+0 
SUMA   =   2.52U01D+0 
TIME   =    0.78988D-1 
ANALYTICAL   =    I.54308D+O 
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Case    18. 


X 


-1 


(i2  +  a2 


sin  at  —  at  cos  at 
2a? 


SL6ROUT 
THIS  IS 
I VPLICI 
I  wplICI 

ZERO    = 
CS1    =    0 
XPTR    = 
XPTI     = 
CXPT    = 
CNUM    = 
CS4    =    C 
CS5   =    0 
CS6    =    C 
CS7   =    C 
CDENO" 
CRESLT 
U   =    ORE 
V    =    DIM 
RETURN 
END 


INE     VALUE     (X,Y ,T,4L,U,V,MYVMM) 

TABLE    ENTRY    NUMBER     18    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 
T    REAL*8     (  A, 8,  C-H,0-Z  ) 
T  C0MPLEX*16     (C) 

0.0+0 

CMPLX(X.Y) 

X*T 

Y*t 

DCMPLX(XPTR,XPTI) 

CDEXP(CXPT  } 

S1*CS1 

CMP LX1AL, ZERO) 

S5*CS5 

S4    +    CS6 

=    CS7*CS7 

=    CNUM/COENOM 

AL(CRESLT) 

AG(CRESLT) 


AL   =    1.25D-1 

T   =    8,E>+0 
AA   =    0.3 
XF    =    -4.11 

YF  =  2-i 
MMMMM  =   22 

LOSC    =   2 

SUM  =   7.70992D+1 

SUMA  =  7.70992D1 

TIME   =    0.51870D-1 

ANALYTICAL  7-70992D+1 
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Case    19. 


X 


-1 


(t»  +  a*)* 


— _  t  sin  at 

—  2o 


SUBRG 
THIS 
IMPLI 
I  fPLl 

ZERO 
CS1    = 
XPT? 
XPTI 
CXPT 
CS3   = 
CNJM 
CSA   = 
CS5    = 
CS6    = 
CS7    = 
CCENO 
CRESL 
U    =    D 
V    =    0 
RETUR 
ENO 


UrINE 

IS    Ti 

CIT    R 

CIT   C 

=    MM 

=     O.D 

DCMP 

=     X*T 

=    Y  *T 

=    OCM 

CDEX 

=    CS1 

CS1* 

DC^P 

CS5* 

CS4 

M     =     r 

=  c 

RE&K 
I*AG( 
N 


VALUE     (X  ,Y  ,T,  AL,UtV,MMHMMJ 

8LE    ENTRY    NUMBER    19    OF     LAPLACE 
EAL*8    ( A, B, D-h,0-ZJ 
OMPLEX*lo     (C) 

MMM     +      1 

+0 
LX(X,  Y) 


■RANSFORMS    BY    SPIEGEL 


PLX(XPT*,XPTI J 

P  <  C  XP  T  ) 

*CS3 

'  SI 

LX(AL .ZERO  ) 

CS5 

♦    CS6 

S7*CS7 

NUM/CDENOM 

CRESLT) 

CRESLT) 


AL    *    1.25D-1 

T    ■    8.D+0 

AA   =    0.3 

XF    =    -4.54 

yF    =   2.2 

MMMMM   =23 

LOSC    =    3 

SUM   =   2.69271D+1 

SUMA   =   2.69278D+1 

TIME  =   0.891S0D-1 

ANALYTICAL  2.69271D+1 
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Case   20 


£ 


"1 


(«*  +  o*)» 


—  sin  at  +  at  cos  e< 

—  2a 


SUB  ROUT 
THIS  IS 
IMPLICI 
IHPLICI 

MM^M^    = 
ZERO    = 
CS1    =    0 
CS2   =    C 
XPTR    = 
XPTI     = 
CXPT    = 
CS3    =    C 
CNUM    = 
CS4    =    C 
CS5    =    D 
CS6   =    Z 
CS7    =    C 
CDENO^ 
CRESLT 
'J    =    DRE 
V    =    DI  M 
RETURN 
END 


INE    VALUE     IXtYtTtALiUtVfMMWMMJ 

TABLE    ENTRY    NUMBER    20    OF     LA°LACE 
T    RE£L*8     (  At  9,  0-t-tO-ZJ 

T    C 3 ^ LEX* 16     (C) 

0.0+0 

cmplx(x,  yj 

S  1*0  S  1 

X*T 
Y*T 

0CMOLX(XPTR  ,XPTI  ) 
DEXD(CXPT) 
CS2*CS3 
S 1  *"  S 1 

CMPLX ( AL.ZERO) 
S5*':S5 
S4    +    CS6 
=    C37*CS7 
=   CMUM/CDEiNOM 
AL(CRESLT) 
AG(C5ESLT) 


TRANSFORMS    BY    SPIEGEL 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =    0.3 
XF    =    -k.$k 

yp  =  2.2 

MMMMM   =   23 

L03C    =  U 

SUM  =   5.52709D+0 

3UMA   =   5.89548D+0 

TIME   =   0.60684D-1 

ANALYTICAL        5-52709D+0 
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Case  21 


X 


-1 


(s2  +  a')* 


__  cos  at  —  ^at  sin  at 


SUBROUTINE    V*,LUE     (  Xt  Y  t  T,  H  ,  U,  V  ,  MMMMfM 

THIS     IS    TA3LE    ENTRY    M'JMBER    21    OF    LAPLiCE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    REAL*8    ( 4, 3 , O-H ,0-Z ) 

IMPLICIT    C1MPLEX*16    <C  ) 

MKMMM    =     M«MMM    +    1 

ZERD    =    O.D+O 

CS1    =    DCMPLX(X,Y) 

CS1    =    DCMDLX(X.Y) 

CS2   =    CS1*CS1*C31 

XPTR    =    x*T 

XPTI     =    Y*t 

CXPT    =    DCMPLX(  XPTR,  XPTI  ) 

CS3    =    CDEXP(CXPT) 

CNUM    =    CS2*CS3 

CS4    =     CS1*CS1 

CS5    =    DCMPLX(ALtZERO) 

CS6    ■    CS5*CS5 

CS7   =    CS4    ♦   CS6 

CDENOM    =    C*7*CS7 

CRESLT    =    :NU,-'/COENOM 

U   =    DREAU  CRESLT) 

V    =    DIMAG(CRESLT) 

RETURN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.3 
XF  =  -4.99 
yP  =  2.3 

MMMMM   =   2k 

LOSC    =   5 

SUM   =    1.19567D-1 

SUMA   =    l.fl94fllD+0 

TIME   =    0.10865D+0 

ANALYTICAL        1.19567D-1 
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Case   22 


X 


-1 


(«*  +  aS)3 


t  cos  at 


SU8*">U 
THIS     I 

impli: 

I"PLIC 
MMMMM 
ZERD    = 
CS1    = 
CS2   = 
XPTR    = 
XPTI     = 
CXPT    = 
CMUM    = 
C^3    = 
CDENOM 
CRESLT 
U    =    DR 
V   =    DI 
RETJRN 
END 


TIME     VALUE     (X,Y,T,iL,U,V,MMMMM) 

S    TABLE    ENTRV     NUMBER    22     OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IT    RE AL*8     (  A,B,D-H,  O-Z) 

IT    Cn"°lEX*lb    (CJ 

a     MVMMM    +     I 

0.0+0 
DCM0LX (X,  Y) 
OCMPLX(AL.ZERO) 

X*T 

Y*T 

DCMDLX(XPTR,XPTI) 

(CCS1*CSU    -     (CS2*CS2)  )*CDEXP(CXPT) 
(  <CS1*CS1)     +    (CS2*CS2)J 

=    C  S3*CS3 

=    CNUM/CDENOM 
EAL(CRESLT) 
MAG(CRESLT) 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =   0.3 
XF   =  -U.54 
yp  =  2.2 

MMMMM   =23 

LOSC    =   U 

SUM  =  U.32242D+0 

SUPIA   =   U.76OOOD+O 

TIME   =    0.67210D-1 

ANALYTICAL        U.322^2D+0 
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Case   23 


£ 


-1 


i 


(S-  -  a2)2 


—  at  cosh  at  —  sinh  a( 

—  la? 


SUBROUTINE    VALUE    (X  ,  Y  ,T  ,  *Lt  J  ,V  ,  MMMMM  ) 

THIS    IS    TABLE    ENTRY    NUMBER    23     OF     LAPLACE    TRANSFORMS 

IMPLICIT    REM_*S     (  A,3,0-H,0-Z) 

I  VPLICIT    C3UPLEXJ<16     (C) 

MMMMM    =     MMMMM    +     1 

ZERO     =    O.O+O 

CS1    =    OCMPLX(X,Y) 

XPT5     =     X*T 

X  FT  I    =    Y  *T 

CXPT    =     DCMPLX(Xt>TR,XPTI  } 

CNUM    =    CDE  XP(CXPT) 

CS<»    =    CSl*CSi 

C  S5    =    OCMPLX(ALtZcRO) 

CS6    =     CS5*CS5 

CS7    =    CS4    -    C36 

CDENOM    =    r  S7*CS7 

CRESLT    =    CNUM/C0EN3M 

U   =   3REAUCRESLT) 

V  =    OIMAG(CRESLT) 

RETJRN 

END 


BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA    =    0.4 

XF    =    -k.kk 
yF   »  2.2 

MMMMM   =23 

LCSC    =   2 

SUM  =   9.^17712+1 

SUMA   =   9.41771D+1 

TIME   =   0.61984D-1 

ANALYTICAL        9.^1771D+1 
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Case  2k. 


£ 


~1 


(«-'  -  a-')* 


t  sinh  of 
2a 


SIB'.OUTINE     VALUE     (  X  ,  Y  ,  T  ,AL,  U  ,  V  ,  mmmmmj 
THIS     IS    TA5LE    ENTRY    MIJM9ES    2<*    n.  r    LAPLACE 
IMPLICIT    REAL*8     (A,3,C-H,>Z) 
IM3L  IC  IT    CT-<PLEX*16    IC  ) 
M^MMM    =     M^MM    ♦    1 

ZERO    =     0.0+0 

CSl    =     DCMPLXU.Y) 

XPTR    =    X*T 

XPTI     =    Y*T 

CX°T    =     OC^LXl  XPTR,  XPTI  ) 

CS3    =    CDEXO(CXPTJ 

CNU*    =    CSl*CSi 

CS4    =    CS1*CS1 

CS5    =    DC"PLX(  AL.ZERO) 

CS6    =    CS5*CS5 

CS7    =    CS4    -    CS6 

COENOM    =   CS  7*CS7 

CRESLT    =    CNUM/CDENHM 

U    =    DREALtC^ESLT i 

V    =    01  MAG(CRESLT) 

RETJRN 

END 


TRANSFORMS    BY    SPIEGEL 


AL   =    1.25D-1 

T   =    8.D+0 

AA   ■    0.4 

XF   =    -4.44 

yF  =  2.2 

MMMMM  =23 

LOSC    =    3 

SUM  =   3.76064D+1 

SUMA   =    3.76151D+1 

TIME   =   0.66924D-1 

ANALYTICAL        3.76064D+1 
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Case   25. 


£ 


-1 


(»-•  -  **p 


sinh  at  +  at  cosh  at 
2a 


SUERGU 
THIS    I 

impli: 

IVPL1C 
MfMMM 
ZERD    = 
CS1    = 
CS2    = 
Xt>TP. 
XPTI     = 
CXPT    = 
CS3   = 
CNUM    = 
CS4   = 
CS5    = 
CS6    = 
CS7    = 
CDENOM 
CRESL7 
U   =    DP 
V    =    01 
RETURN 
END 


TIME    VALUE     (X,Y,T, aL,U,V,VMMHM) 
S    TABLE    ENTRY     NUMBER    25    OF    LAPL4C 
IT    PEAL*3     (  A,8,D-H,3-Z) 
IT    COMPLEX* 16     I C  J 
■     MMMMM    +     1 

0.0+0 
OCMPLX(X,YJ 
CS1*CS1 

=      X*T 

Y*T 

DCHPLX(XPTR,XP7I) 
CDEXD(CXPT) 

C<:2*CS3 
CS1*CSI 

DCMPLXULtZEROJ 
CS5*CS5 
CS4  -■ CS6 

=    CS7*CS7 

=   CNUM/COENOM 
EAL(CRESLT) 
MAG(CRESLT) 


E    TRANSFORMS    BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    O.U 
XF    =   -4.89 

yF  =2.3 

MMMMM  =   2h 

LOSC    =  k 

SUM  =   1.08731D+1 

SUMA   =    1.20870D+1 

TIME   =   0.73710D-1 

ANALYTICAL   =    1.08731D+1 
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Case   26. 


X 


-1 


(«2  _  a2 


cosh  at  +  latsinhat 


SUBR? 

THIS 

IMPLI 

IM°LI 

MMMMM 

ZERO 

CS1   = 

CS2    = 

XP~R 

XPTI 

CXPT 

CS3    = 

CMUM 

CS4    = 

CS5   = 

CS6    = 

CS7    = 

CDEM3 

CRESL 

U   =    D 

V   =   3 

RETJR 

END 


UTINE     VALUE     (  X  ,  Y  ,  T  ,4  L,  'J ,  V  ,  MMMM") 
IS    TARLE    E\7RY    NUMBER     26    IF     LAPLACE 
CIT    REAL*8     { AfB,C-H, 0-Z) 
(C) 


TRANSFORMS    BY    SPIEGEL 


CIT   COMPLEX*  16 

=       MM-JMM       +       I 

=    0.0+0 

OCMPLXC  X,Y) 

CSI*CS1*CS1 
=    X*T 
=     Y*T 
=     DC^PLX  (XPTR,XPTI  ) 

CDEXP(CXPT) 
=    CS2*CS3 

CS1*C31 

DCMPLX(AL,ZERO) 

CS5*CS5 

CS*    -    C36 
M    =   C  <-7*CS7 
t    =    CNU^/CDtNOM 
REALICRESLT) 
IMAG(CRESLT) 
N 


AL   =    1.2  5D-1 
T   =   8.EH-0 
AA   =   0.4 
XF  =  -4.89 

yF  =  2.3 

MMMMM   =24 

LCSC    =   5 

SUM  =   2.13068D+0 

SUMA   =   4.30908D+0 

time  =  0.19490D+0 

ANALYTICAL        2.13068D+0 
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Case   27 


X 


-1 


t  cosh  a  t 


SU8R3 

THIS 

I  MPLI 

IMPLI 

MMMMM 

ZERO 

CS1    = 

CS2    = 

XPTR 

X°TI 

CXPT 

CS3  = 
CDENO 
CRESL 
U  =  D 
V  =  3 
RETUR 
END 


UTINE 
IS  TA 
CIT  P 
CIT  C 
=    yv 

=    0.0 
DCMO 

OCMP 
=  X*T 
=  Y*T 
=  CCM 
=     <IC 

(  (CS 
M  =  C 
T  ■  C 
R  E  AL  ( 
IMAG( 


VALUE     ( X, Y,T, AL,U,V,MyMMM) 
8LE    E4TRY    NUMSER    27    OF    LA-.°LACE    TRANSFORMS    3Y    SPIEGEL 
EAL*8     ( A, 3tD-h,3-Z) 
0'1PLEX*lb    (C) 
MMM    +    1 
+  0 

LX(  X,  Y) 
LXULiZERT} 


DLX(XPTR,XPTI  ) 

S1*CSU     +    (CS2*CS2))*CDEXP(CXPTJ 

1*CSU    -    (CS2*CS2>) 

S3*CS3 

NUM/C9EN0H 

f^ESLT) 

CRESLT) 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.4 
XF    =   -4.89 

yp  ■  2«3 
mmmmm  =  24 

LOSC    =   4 

SUM   =    1.23446D+1 

SUMA  =   1.3I*023D+1 

TIME   =   0.65078D-1 
ANALYTICAL        1.23446D+1 
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Case   28. 


X 


-1 


(«S  +  a-')3 


(.1  —  a2t-)  sin  at  —  3«<  cos  at 
ia- 


SUBROUTINE     VALUE     (Xt  Y»T,  4LtUf  V,  MMVMVJ 

THIS    IS    ""ISLE    E^TRY    NUMBER    28    IF     LAPLACE 
IMPLICIT    PEAL*3     ( A, 3, C-h,0-Z) 
IMPLICIT   C-implEX*16    (C) 

MVM^M      =       MMMMM      +       I 

ZERO    =    O.O-t-0 

CS1    =    DCMPLX(XtY) 

XPTR     =     X*T 

XPTI     =    Y*T 

CXPT    =    DCMPLXt  XPTR,XOTI  ) 

CNUM   =    CDEXP(CXPT) 

CS4   =    CS1*CS1 

Z  S5    =    DCMPLX(ALfZERO) 

CS6    =    CS5*CS5 

C  S7   =    CS4   +    CS6 

CDE"nM  =  rs7*r.s7*,:s7 

CRESLT    =   CNUM/C3cN0M 
U    =    JRE1LI  CRESLT) 
V    =    DIMAG(CRESLT) 
RETURN 
END 


TRANSFORMS    BY    SPIEGEL 


AL   =   1.25D-1 

T  =   8.D+0 

AA   =    0.4 

XF  =  -3.6 

yF   =   2.0 

MMMMM   =21 

LOSC    =■  2 

SUM   =   2.54096D+2 

SUMA   =   2.89148D+2 

TIME   =   0.5041UD-1 

ANALYTICAL   =   2.54096D+2 


120 


Case   29 


X 


-1 


(*-  +  «-•)■•' 


t  sin  at   —  nt2  cos  a t 
8uJ" 


SUBRCU 
THIS  I 
I VPLIC 
I MPLIC 

ZERO    = 
CS1    = 
XPTR    = 
XPTI     = 
CXPT    = 
CS3    = 
CMJM    = 
CS4   = 
CS5    = 
CS6    = 
CS7   = 
COENOM 
CRESLT 
U    =    0° 
V    =    01 
RETURN 
ENO 


TINE     VALUE     (X , Y  ,  T  ,  ALtU , V , «mmv v ) 

S    TABLE    ENTRY    NUMBER    29    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 
IT    RE£L*3     (  At8,D-F,">-Z  ) 
IT    COMPLEX*  16     tC) 
=     MM^MM     +     I 

0.0+0 
DC^PLX(  X,Y) 

X*T 

Y*t 

DC"PLX(XP7R,XPTI  ) 
COEX  °(CX=T  ) 

CS1*CS3 
CS1*CS1 

OCMDLX  (oL,  ZERO  ) 
CS5*CS5 
CS4   -t-    CS6 

=    CS7*CS7*CS7 

=   CNUM/COE  NOM 
EAL(C'ESLT) 
MAG(CRESLTJ 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.3 
XF    =    -4.11 

yF  ■  2.1 

MMMMM   =   22 

LCSC    =   3 

SUM   =    1.5U198D+2 

SUMA  =  1.66758D+2 

TIME    =    0.56U72D-1 
ANALYTICAL   =    1.5^198D+2 
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Case  30. 


X 


-1 


(«-  +  a-)3 


""""  (1  +  g-l-)  sin  at  —  at  ens  at 


SUBROUTINE    VALUE     (  X,  Y,  T  ,ALt  U,  V,  mmmmm) 

THIS     IS    TABLE    ENTRY    NLMQE'    33    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

I  MPLi:  IT    REAL*8     (  A  ,  8  ,  D-h  ,0-Z) 

IMPLICIT    CTMPLEX*I6    <  C  ) 

MMMMM    =     ^MMMM    +    1 

ZERO    =    O.D  +  O 

CS1    =     DC*PLX(X,Y) 

CS2    =    CS1*CS1 

XPT*    =     X*T 

XOTJ     =     Y*t 

CXPT    =    DC^PLX(XPTR,XPTI) 

CS3   =    CDEXP(CXPT) 

CNUM    =    CS2*CS3 

CS4   =    CSI*CS1 

CS5    =    0C^PLX(AL  tZERO) 

CS6    =    CS5*CS5 

C  S7    =    CS4   +    CS6 

CDEMDM    =    CS7*CS7*CS7 

CRESLT    =    CNUM/CDEMO'1 

U    =    DREAL(CRESLT) 

V    =     DIM AG(CRESLT) 

RETURN 

ENO 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =    0.3 
XF   =   -4.11 

yF  =  2.1 

mmmmm  =  22 

LOSC    =   2 

SUM   =   7.312S9D+1 
SUMA   =   7.31289D+1 
TIME   =   0.57122D-1 
ANALYTICAL   =    7.31239D+1 
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Case   31 • 


£ 


-1 


(«-  +  a-)1 


31  Bin  nt    t-   at-  ro*  at 
8a 


SUP."  3 
ThIS 
IMPLI 
IMOLI 

MMMMM 
ZERO 
CS1    = 
CS2    = 
X°T3 
XPrI 
CXDT 
CS3    = 
CNUM 
CS4    = 
CS5    = 
CS6    = 
CS7    = 
CDENO 
CRESL 
U  =    0 
V    =    0 
RETUR 
END 


IN  5  V 
TARL 
T  ^  E\ 
T  COM 
MMMM 
O.D  +  0 
CMOLX 
S1*CS 
X*T 
Y*T 
DCMPL 
DEX^! 
CS2*C 
S1*CS 
CM:>LX 
S5*CS 

=    C<^  7 

=  GNU 
A  L  (  C  R 
AG(CR 


ALJE     (X,Y,T,4Lf Ut Vt^MM^M) 

E    ENTRY     NUMBER    31    IF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

L*8    (  A, 3  ,0-H,  3-D 

°L  EX*  16    (  C  J 

'■<   +    1 

(X,  YJ 
l#CSl 


X(XPTR,XPTI ) 

CXPT) 

S3 

1 

< AL.ZERO) 

5 

CS6 

*CS7*CS7 

M/CutNOM 

ESLTJ 

ESLT) 


AL   =    1.25D-1 

T  =   8.D+0 

AA   =    0.3 

XF  =  -b.si* 

yF    =   2.2 

HMMMM   =23 

LOSC    =    3 

SUM   =   2.U5177D+1 

SUMA  =  2.U<l84D+l 

Tire  =   0.6055UD-1 

ANALYTICAL   =   2.^5177^+1 
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Case   32. 


£ 


-1 


(j2  +  a- 


— ■  (3  —  a-t-)  sin  at  +  Sat  cos  at 
8a 


SLSOUTIN5     VALUE     (  X  ,  V  ,7  ,  AL  t  U  t  V  f  «MH«m) 

THIS     IS    TABLE    EnIT'.Y    NUM9E*     32    DF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    PEAL*3     ( A, 3, C-h,0-Z) 

IMPLICI  T   COMPLEX*  lo    (C  ) 

PMMMM     =     MMMMM     +     1 

ZERO    =    0.0+0 

CS1    =    OCMPLX(X,Y) 

CS2    =    CS1*CS1*CS1*CS1 

XPTR    =    X*T 

XPTI     =    Y*T 

CXPT    =     DC^°LX  (XPTR,  XPT  I) 

CS3    =    CDEXP(CXPT) 

CNLM    =     CS2*TS3 

CS4    =    CS1*CS1 

CS5    =    DC*DLX<  ALtZERO) 

CS6   =    CS5*TS5 

CS7    =    CS4    +    C36 

:denom  =  :s7*cs7*cs7 

CRESLT    =    CN'JM/COENOM 
U    =    3REAL(0RESLT) 
V   =    DIMAGCRESLT) 
RETJRN 
END 


AL   =    1.25D-1 

T   =   8.D+0 

AA   =   0.3 

XF    =   -^.5^ 

yF    =   2.2 

MMMMM   =   23 

LOSC    =   i* 

SUM   =   ^.38^^5D+0 

SUMA   =   iv.830^8D+0 

TIME   =   6.7548D-2 

ANALYTICAL  =  U.38^5D+0 
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Case   33 


£ 


-1 


(« s  +  a-)-' 


(8  —  a-i-')  cos  a f  —  Tat  sin  at 
8 


SUBROUTINE    VALUE    (X  ,Y  ,T  ,  AL,  J  ,V  ,  «M«mm  ) 
THIS    IS    T6=SLE    ENTRY    NUMBER    33    OF     LAPLACE 
IMPLICIT    REAL*8     (A,8,D-h,>Z) 
I  VPLICIT    COMPLEX*  16     (C) 

M'WMM     =     MMMMM     +      I 

ZERO    =    O.O+O 

CS1    =    OCMPLX(X,Y) 

cs2  =  ;si*':si*csi*csi*csi 

XPTR    =    X*T 

XPTI     =     Y*T 

CXPT    =     DC^fLXt  XPTR,  XPTI  ) 

CS3    =    CDEX°(CXP-) 

CNUM    =    CS2*CS3 

CS*    =     CS1*CS1 

CS5    =    DCMPLX(AL, ZERT) 

CS6    =    CS5*CS5 

CS7    =    CS4    +   CS6 

cdenom  =  cs7*cs7*c?7 
creslt  =  :num/cdenom 

u  =  dreal(c0e3lt) 
v  =  oimag(creslt) 

RETJRM 
ENO 


TRANSFORMS  BY  SPIEGEL 


AL  =  1.25D-1 

T  =  8.D+0 

AA  =  0.3 

Xj,  =  -4.99. 

yF  =  2.3 

MMMMM  ■  2& 

LOSC  =  5 

SUM  =  -2.63523D-1 

SUMA  =  1.8l9^3D+0 

TIME  =  6.7MH.1+D-2 

ANALYTICAL   =    -2.63523D-1 
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Case  3k 


£ 


-1 


t-  sin  at 
2a 


SUBROU 
THIS  I 
IfPLIC 
I  fPLIC 

ZEP.O   = 
CS1    = 
CS2    = 
XPTR    = 
XPTI    = 
CXPT    = 
CNUM    = 
CS3   = 
CDENOM 
CRESLT 
U    =    D= 
V    =    01 
RETURN 
END 


TINE     VALUE     (X tY ,T , AL,U , 7 , "mmmv ) 

S    T43LE    ENTRY    NUMBER    34    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IT    REAL*8     (  A, 3,  D-h,0-ZJ 

I  T   CQmdlEX*16     (  C) 

=      MMMMM     +     I 

O.D+0 
OCMPLXJ  X,  Y) 
0CMOLX(ALi2ER0 J 

X*T 

Y*T 

OC«PLX(XPTR,XPTI ) 

( (3.0+0* (CS1*CS1)1-(CS2*CS2) )*COEXP(CXPT  ) 
(CS1*CSU    •■    (CS2*CS2) 

=    CS3*CS3*CS3 

=   CNIJM/CDENOM 
EAL(CP.ESLT) 
MAG(CRESLTJ 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =   0.^ 
xF   =   -k.kk 

yF  =  2.2 

MMMMM   ■   23 

LOSC    =   2 

SUM  =  2.15417D+2 

SUr/IA   =   2.15417D+2 

TIME  =  0.6783UD-1 

ANALYTICAL   =   2.I5U7D+2 
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Case  35. 


<£ 


-1 


3n-s 


mm „  £*2  cos  at 


SUBROUTINE    V&LUE    (  X  ,  Y  ,T,A  L,U,V  ,  MMMMV) 

THIS     IS    TABLE    ENTRY    NUMBER    35    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICI  T   RE&L*3    (  A  ,8  ,0-H  ,0-Z) 

IMPLICIT    CTMPLEX*16     (C ) 

MMMMM    =     ^mmmm    ♦    1 

ZE^Q    =     0.3-1-0 

CSl    =     DC**°LX(X,Y) 

CS2    =    OCMPLX(AL.ZERT) 

XPTR    =    X*T 

XPTI     =     Y*T 

CXPT    =    DCMPLX  (XPTR,X»TI  ) 

CNUM    =     (<C  31*CSl*CSl)-(  2.D+0*CS1*CS2*CS2)  )*CDEXP(CXPT) 

CS3    =     <CS1*CS1)     +    <CS2*CS2) 

COENOM    =    CS3*CSi*CS3 

CRESLT    =    CNUM/CDEN04 

U    =    DREAL(CPESLT) 

V   =    OIMAG(CRESLT) 

RETURN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =   0.3 

xF  =  -4.5U 

yF  =  2.2 

MMMMM   =   23 

LOSC    =   3 

SUM   =   1.72897D+1 

SUMA   =   1.72904D+1 

TIME   =   8.447/+D-2 

ANALYTICAL   =    1.72S97D+1 
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Case  36 


£ 


-1 


s*  —  6a-V  ■*-  a4 
(a*  +  a')* 


SUBROUTINE     VALUE     (  X  ,  Y  ,  T  ,  41, U,  V  ,  MMMMM) 

THIS     IS    TARLE    E'JTaY    NUMBER     36    ">F     LAPLACE     TRANSFORMS    BY    SPIEGEL 

IMPLICIT    REAL*3     ( A, 8, C-H,3-Z) 

IMPLICIT   COMPLEX*  16    (C) 

MMMMM    =     MMM'^m    +    1 

ZERO    =    O.D+O 

CS1    =    DCMDLX(X,Y) 

CS2    =    DC'-IPLX  (  AL,ZE'0) 

CS3    =    CS1*CS1*C>1*C~1 

CS4   =     6.D+0*T5^*:  S2*CS1*CS1 

CS5    =    CS2*CS2*CS2*CS2 

XPTR    =     X*T 

XPTI    =     Y*T 

CXPT    =    DCMPLX(XPTR,XPTIJ 

CNUM    =     (CS3-CS4-+:  S5)*CDEXP(CXPTJ 

CS6    =     <(CS1*CS1)     +    (CS2*CS2)) 

CDENOM    =    CS6*CSo*Cc6*CS6 

CRESLT    =    CNUM/COENOM 

U    =    DPEAL(CRESLT) 

V    =    DI MAG(CRESLT) 

RETURN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =0.^ 
xF    =   -b.kU. 
yF  =  2.2 

MMMMM   =   23 

LOSC    =   2 

SUM   =   1*.61058D+1 

SUMA   =  ^.6l058D+l 

TIME   =    0.73^2UD-1 

ANALYTICAL   =   ^.6l058D+l 


128 


Case   37 


e? 


"I 


g3  —  a2s 
(»-  +  a-)4 


t3  sin  at 
24a 


SUBROUTINE     V&LUE     (X,Y,TfAL, 

U, 

V  ,MMM<4M) 

THIS     IS    TABLE    E-JTRY    NUMBER 

37 

OF 

LAPLACE 

TRANSFORMS 

IMPLICIT    REAL*3    (A,3,D-Ht0- 

■Li 

IWPLICIT    r.H  «t>LEX*16    (C  ) 

MMHMM    =    MMMMM    +    L 

ZERO    =     0.0+0 

CS1    =     DCMPLX(X,Y) 

CS2    =    DCMPLX(ALtZERD) 

XPTR    =    x*T 

XPTI    =     Y*T 

CXPT    =    DCM  =  LX  (XPT  R,XPTI  ) 

CNUM    =     <(CSI*CS1*CS1)     -    (CS2* 

CS2* 

C31)  )*CDE 

CS3    =     <(CS1*CS1)     +    (CS2*CS2)) 

CDENO*    =    CS3*CS3*CS3*CS3 

CRESLT    =   CNUM/CJEN3M 

U    =    OREAKCESLTJ 

V   =    DIMAG(CR-ESLT) 

RETURN 

END 

BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =   O.U 
xF  =  -4.01 

yF   -  2.1 

MMMMM   =   22 

LOSC    =   3 

SUM  =  1.43611D+2 

SUMA   =l.U4005D+2 

TIME  =  0.68692D-1 

ANALYTICAL  =  1.^36llD+2 
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Case  38. 


J?'1 


(,2  _  a2)3 


(3  +  aU*)  sinh  at  -  Sat  cosh  at 

8a* 


SU8R0U 
THIS  I 
IMPLIC 

IfPLIC 
MMMMM 
ZERO    = 
CS1    = 
XPT*    = 
XOTI    = 
CXPT    = 
CNUM    = 
CS4    = 
CS5    = 
CS6   = 
CS7    = 
CDENOM 
C&ESLT 
U    =    OR 
V    =    01 
RE^JRN 
ENO 


TINE    VALUE     (X,Y,T,  iL,U,  V.^MMMM) 

S    TABLt    t^TRY    NUMBER    38    OF    LAPLACE    TRANSFORMS 

IT    ?  EAL-8    ( A,R,0-H,0-Z ) 

IT    C3'1=>LEX*lo     (  CJ 

=     MMMMM    +    1 

0.0  +  0 
DC'^LXU,  YJ 

Y*T 

CC^PLX  (X?TR,X»T  I) 

CDE  XO(CXPT) 
CS1*C>  1 

OCMPLX (AL.ZEROJ 
CS5*CS5 
CS4    -    CS6 

=    C*  7*CS7*CS7 

=    CNUM/CDENOM 
EAL(C=ESLT J 
MAG(CRESLT) 


BY    SPIEGEL 


AL   =    1.25P-1 

T   =   8.D+0 

AA   =    O.U 

Xp   =   -4.01 

yF   =   2.1 

MMMMM  =22 

LOSC    =   2 

SUM   =   2.93122D+2 

SUMA   =   3.6?i4.?lD-t-2 

TIME   =    O.U833^D-l 

ANALYTICAL   =   2.93122D+2 
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Case   39, 


X 


-l 


(«2  -  c2)-1 


at-  cosh  at  —   tsinhat 
8a3 


subou 

THIS    I 

i  mo li  : 

I VPLIC 
MfMMM 
ZERO    = 
CS1    = 
XPTR    = 
XPTI    = 
CXPT    = 
C  S3    = 
C'^JUM    = 
CS4    = 
CS5   = 
CS6    = 
CS7    = 
CDEN3M 
CCESL'- 
U   =    OR 
V    =    31 
RETURN 
END 


TINE    VALUE     (XtY,T,AL»l.J,V,WMMMM) 

S    "ABLE    ENTRY     NUMBER    39    C F    LAPUCE   TRANSFORMS    BY 

IT   RE*L*8     ( A,8,Q-HfO-Z) 

jr    C^'^LEXMo    (CJ 

C.D+O 
DCM°LX(X,Y) 

X*T 

Y*T 

DCMDLX  (X°TR,XOTi  ) 
CDEXP(CXPT) 

CS1*CS3 
CS1*CS1 

DCMPLX(ALt  2EROJ 
CS5*C^5 
CS<t    -    C36 

=   CS7-CS7*CS7 

=    CNUM/COENGM 
EAL(CRESlT) 
MAG(CRESLT) 


SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.4 
XF    =   -4.01 

yp    =   2«1 
MMMMM  =   22 

LCSG    =    3 

SUM  =  1.8835UD+2 

SUMA  =  1.89195D+2 

TIME  =  0.4875OD-1 

ANALYTICAL  =  1.8835^0+2 
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Case   kO, 


£ 


(^  -  a-')3 


at  cush  at  +   la-t-  -  1)  sinh  at 
8a3 


SIBRDUTINE     VALUE     (  X  ,Y  ,T  ,AL,U, V  ,  VMMvm) 

THIS     IS    TiPLE    EWTRY    NUMBER    4-0    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    REAL*3     (A»8tO-h,0-ZJ 

IMPLI2I  T   COMPLEX*  16     {  C) 

M^mmm    =    mm^mm    +    l 

ZERO    =    0.0+0 

CS1    =    DCMPLX(  X,  Y) 

CS2    =    CS1*CS1 

XPTR    =    X*T 

XPTI     =     Y*T 

CXPT    =     DCMOLX(XPTRtXt>TI  ) 

CS3    =    CDEXP(CXPT) 

CNUM    =    CS2*<"S3 

CS*    =    CS1*CS1 

CS5   =    0CMDLX(AL,ZER0) 

CS6    =     CS5*CS5 

CS7    =    CS4    -    CS6 

CDE^DM    =   C<^7*CS7*CS7 

CPESLT    =    CNIJM/CDEN'IM 

U    =    DREAL(CRESLT) 

V    =    DIMAG(CSESLT) 

RETURN 

END 


AL  =    1.2  5D-1 
T  =    8.D+0 
AA   =    0.4 
Xp    =   -4.44 

yF  =  2.2 

MMMMM   =23 

LOSC    =   2 

SUM  =   9.87572D+1 

SUMA   =   9.87572D+1 

time  =  0.5252D-1 

ANALYTICAL   =   9.87572D+1 
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Case  k] 


X 


-1 


(S--a-)* 


— ""  31  sinh  at  +   at2  cosh  at 
8a 


THIS  IS 
IMPLlf  I 
IVPLICI 

MfMMM    = 
ZER3    = 
CS1   =    o 
CS2   =    C 
X°TR    = 
XPTI    = 
OPT    = 
CS3   ■     C 
CNUM    = 
CS4  =    C 
CS5    =     D 
CS6    =    C 
CS7    =    C 
CDENOM 
CRESL"1- 
U    =    ZIRE 
V    =    DIM 
RETURN 
END 


TABLE    ENTRY     NUMBER    41    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 
T   REAL-8    (  Ai3t  D-Hf-3-Z) 
T    COMPLEX* 16    (CJ 

M.MMMM    *■     1 
O.D+0 

CMPLX (X, Y) 
S1*C.S1*CS1 
X*T 
Y*T 

DCM«>LX(XPTR,XPTIJ 
DEXP(CXPT) 
C^2*C*3 
S1*C  SI 

CMDLXUL  ,ZER0  ) 
S5*C35 
S4   -    CS6 
=    CS7*CS7*CS7 
=   CNUM/COENOM 
AL(CRESlT) 
AG(CRESLT) 


AL   =    1.25D-1 

T   =   8.D+0 

AA   =   O.U 

Xp    =   -4.W 

Iff    =   2.2 

MMMMM   =   23 

LOSC    =   3 

GUM   =  4..05U95D+1 

SUMA   =  U.05582D+1 

TIMS  =   0.53508D-1 

ANALYTICAL   =   k.05U95D+l 
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Case  hi 


£'x 


(,2_a-)i 


(3  +  a;<-)  sinh  at  +  hat  cosh  at 
8a 


SUBROUTINE     VALUE     (  X  ,Y  ,T  ,  AL,U  ,  V  ,  MMWmm) 

ThIS     IS    TiBLE    EMTRV    NUMBER    42    DF    LAPLACE 

IMPLICIT    REAL*9     (  A,  B  ,  0-n,  1-1 ) 

IMPLICIT  -OMPLEX*  L6    (  C) 

MMMMM     =     MMMMM    +     l 

ZERO    =    0.0+3 

CS1    =    DCMPLXI  X,Y) 

CS2    =    CSl*CSl*CSi*CSl 

XP^R    =    x*" 

XPTI     =     Y*T 

CXPT    =     DCMOLX  (XOTR,  XPTI  J 

CS3    =    CDEXP(CXPT) 

CNUM    =    CS2*r.S3 

CS4    =     C«U*CS1 

CSS   =    DCMPLXIALtZEROJ 

CS6    =     CS5*'"S5 

CS7    =    CS4-    -    CS6 

CDENOM    =   CS7*CS7*CS7 

CRESLT     =    CNUM/CDENOM 

U   =    DREAL(CRESLT) 

V    =    0  I  MAG(  CRESLT) 

RETURN 

END 


TRANSFORMS    BY    SPIEGEL 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =    Q.k 

Xj,    =    -4.89 
Vp    -   2.3 

MMMMM   =   2k 

LOSC    =  4 

SUM   =    1.2M62D+1 

SUMA   =   1.34824D+1 

TIME   =   0.58240D-1 

ANALYTICAL   =    1.24162D+1 
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Case  43. 


£ 


-1 


<•*  -  a^P 


(8  +  a"t-)  cosh  at   +   1atz\i\hat 
8 


SUBROUT 
THIS  IS 
INPLICI 
I MDLICI 

MMMMM     = 

ZERO    = 
CS1    =    0 
C$2  *    C 
XPTR    = 
XPTI     = 
CXPT    = 
CS3    =    C 
CNLM    = 
CS4    =    C 
CS5    =    D 
CS6   =    C 
CS7    =    C 
CDENOM 
CR6SLT 
i|    =    ORE 
V    =    DIM 
RETURN 
END 


I  ME    VALUE     (XfY»T,ALiUiVt  WMM) 

TABLE    EMTRY    NIJM9ER    43    OF    LAPLACE 
T    <?EAL*8     <  A,  S,  C-l-,0-Z) 
T    COMt>LEX*16     (CJ 

O.D+0 

CMPLX(X,  YJ 

Sl*CSi*CSl*CSi*CSl 

X*T 

Y*T 

CC^OLX(X°TR  ,XPTI  ) 

DEXO(CXPT  J 

CS2*CS3 

S 1  *r  S 1 

CMP  IX  (AL.ZERO) 

S5*CS5 

S4    -    CS6 

=    CS7*CS7*CS7 

=   CNUM/C3EN0M 

ALCCRESL7J 

AC(CRESLT) 


TRANSFORMS    8Y    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.4 

Xp    =    -4.89 

yF  =  2.3 

MMMMM   =   2U 

LCSC    =    5 

SUM  =   2.76427D+0 
SUMA  =   U.7326ID+O 
TIME   =    0.62920D-1 
ANALYTICAL   =   2.?6U2?D+0 
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Case  hk 


X 


-1 


3s-  +  n- 
(u-'-a-)1 


t2  sinh  at 
2a 


SUBROUT 
THIS  IS 
IMPLICI 

INPLICI 
•4MMMM     = 
ZE«0    = 
C51    =    D 
CS2    =    0 
XPTR    = 
XPTI     = 
CXPT    = 
CNJM 
CS3 
CDEMOM 
CR£SLT 
U   =    ORE 
V    =    DIM 
RETURN 
END 


=    ( 


I  ME    VAL'JE    (XfYfT,ALtU»VtMMMMM) 

TABLE    ENT^Y    NUMBER    44    OF     L4°LACE 
t    c  E*L*3    (  A, 6  ,0-H,  >Z  ) 
T    CD^PLcX* 16     ( C) 

MMMMM    +     1 

CM0LX(X,Y) 

CMPLX(AL,7_ER0) 

X*' 

Y*T 

0C'-1OLX(XoTR,XOTI  ) 

(  <3.D+0*(CS1*CS  1)  )+(CS2*CS2)  )*COEXP(CXPT) 

CS1*CS1)    -    (CS2*CS2) 

=    CS3*CS3*CS3 

=    CN'JW/cqsntm 

AL(CRESLT) 
AG(CRESLT) 


RANSFORMS    BY    SPIEGEL 


AL   =    1.25D-1 

T  =   8.D+0 

AA   =   O.U 

Xp    =    -4.M* 

yF   =  2.2 

MMMMM   =23 

LOSC    =   2 

SUM  =  3.00852D+2 

SUMA   =    3.00S52D+2 

TIMS   =    0.550^2D-1 

ANALYTICAL  =    3.00852D+2 
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Case  45 


X'1 


3n-t 


__  |f2  cosh  ot 


S18R0UTINE  VALUE  (X  ,Y  ,T , ALtUt V , mmvm^ j 

THIS  IS  T5.BLE  ENTRY  NUMBER  45  OF  LAPLACE  TRANSFORMS  BY  SPIEGEL 

IMPLICIT    se;l*3    ( A, 8, C-h,1-Z) 

I  KPLICIT    COMPLEX*  lo     (CJ 

MMM"4M     =     MMMMM    +     1 

ZERO    =    O.D+0 

CS1    =    DCMDLX(  X,  YJ 

CS2    =    0CM5LX(4l,ZEHJ> 

XPTR    =    X*T 

XPTI     =    Y*T 

CX»T    =     DCMPLX(XPTR, XPTI J 

CMM    =     (  (CSl*C31*Cil  )+(3.D+0*CSl*CS2*CS2JJ*C0EXP(CXPT) 

CS3   =     (CS1*CS1J     -     (CS2*CS2) 

COENOM    =    CS3*CS3*:S3 

CRESLT    =    CMUM/COcNOM 

U  =    DREAL(CRESLT) 

V    =    DIMAG(CRESLT J 

RETURN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA  =   0.4 

Xp   =   -k.kk 

yF  =  2.2 

MMMMM   =23 

L05C    ■   3 

SUM  =  4.93786D+1 

SUMA   =   4.93892D+1 

TIME   =   0.57772D-1 

ANALYTICAL   =  4.93786D+1 
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Case  46. 


£ 


-1 


s*  4-  Cn-s-  +  a* 


~~ "  |f3  cosh  a< 


SUBR3UT 
Tril  S  IS 
IMPL EC  I 

IMPLICI 
MMMMM  = 
ZE'3    = 


CS1 
CS2 
CS3 
CS4 

:s5 

XPTR    = 
XPTI     = 
CXPT    = 
CNJM    = 
CS6   =    ( 
COEMDM 
CRESLT 
U   =    OPE 
V    =    DIM 
RETURN 
END 


INE    VAl'JE     (XtY,T,ALt'JfV,MMMMM) 

T49LE    EMTRY    NUMBER    4t>    OF    LAPLACE    TRANSFORMS     8Y    SPIEGEL 
T   o  E4L*9    ( At3t0-H,l-Z) 

T   C0'^PLEX*16    (C) 

MM^MM      +■       J. 

0.0+0 

CMPLX(X,Y) 

CMP  IXUL, ZERO) 

Sl*CSl*Cal*CS  1 

,0+0*CS2*CS2*CSl*CSl 

S2*CS2*CS2*CS2 

X*T 

Y*T 

DCMPLX(XPTR,  XPTI  ) 
(C33+CS4  +  CS5 )*CDEXP(CXPT  ) 
(CS1*CSU    -    (CS2*C52)) 

=    CS6*CS6*CS6*CS6 
=   CN'JM/C06N0M 
AUCRESLT) 
AG(CRESLT) 


AL   =    1.25D-X 

T   =   8.D+0 

AA   =    0.4 

xF  =   -4.44 

yF  =  2.2 

MMMMM  =   23 

LOSC    =   2 

SUM  =  1.31676D+2 

SUMA  =  1.31676D+2 

TIME  =  O.65676D-I 

ANALYTICAL  =  1.3l6?6D+2 


138 


Case  k7. 


J?-1 


n1   f  g'.t 


t1  sinh  nf 
21a 


SUBROUTINE     VALUE     (  X  ,  Y  ,T  ,  AL,'J,  V  ,  MMVMM) 

THIS     IS    "ABLE    EMTRY    NUM3ER    47    OF    LAPLACE    TRANSFORMS    8Y    SPIEGEL 

IMPLICIT    REAL*9     ( A, 3 , 0-Hf D-ZJ 

IMPLICIT    COMPLEX*  16    (  C) 

M^MMM    =     M^M"^M    +    l 

ZERO    =     O.D+0 

C  SI    =    OCMPLX(  X,Y) 

CS2    =    DCMPLX(AL, ZER3) 

XPTR    =    x*"1" 

XPTI     =    Y*T 

CXPT    =     DC^PLX  (XPTR, XPTI) 

CNUM    =     ({CS1*CS1*CS1J     +    <CS2*CS2*CS1  )  )*CDEXP(CXPT) 

CS3   =    ((CS1*CS1)    -    (CS2*CS2)J 

CDEMOM    =    CS3*CS3*CS3*CS3 

CRESLT    =    :NUM/CD£NOM 

U   =    DREAL(C3ESLT) 

V    =    DI^IAG(CRESLT) 

RETURN 

END 


AL  =    1.25D-1 
T  =  8.D+0 
AA  =   0.^ 
XF  =   -fc.01 

yF  =  2.1 

MMMMM   =22 

LOSC   =   3 

SUM  =  2.00568D+2 

SUMA  =  2.01796B+2 

TIME  =  0.55822D-1 

ANALYTICAL  =  2.00568D+2 
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Case   1*8. 


X 


-1 


rmri   \    _  \f\nt  \ftot 

—— r  i  \'3  sin  — 3 cos  — r, —   + 

.la-   \  2. 


«l 


SUBROUTINE    VALUE     (X  ,Y  ,T  ,  iL  ,U  ,V  ,  mmmmm  j 

THIS    IS    TA.RLE    ENTRY    NUMBER    -+8    DF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT     REALMS     (  A,B,  C-h,0-Z  ) 

I MPLICIT    C0"PLEX*16     (C) 

MMMMM     =     MM«"M     *      1 

ZERO    =    O.D+0 

CS1    =    OCMPLX(X,Yi 

XPTR    =     X*T 

XP~I     =    Y*T 

CXPT    =    DCMPLX(XPTRtXPTIJ 

CNUM    =     CDEXP(CXPT ) 

CS4    =    CS1*CS1*CS1 

CS5    =    OCMPLX(AL,ZERO) 

CS6    =    CS5*CS5*CS5 

COENOM    =   CS4    +    CSo 

CRESLT    =    :  NUM/CDENOM 

U    =    0REiL(C°  ESL7  ) 

V    =    OIMAG(CRESLT) 

RETURN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    O.U 
Xp   =   .4. 44 

YF    "   2-2 

MMMMM   =   23 

LOSC    =    3 

SUM   =    3.1U683D+1 

SUMA   =    3.1^767D+1 

TIME  =    0.U8958D-1 
ANALYTICAL   =   3.14683D+2 
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Case   k3, 


£ 


-1 


c""  f 


VZat 


,3  +  a' 


SU8R0UTINE    VALUE     (X  ,  Y  ,  T  ,  AL,U  ,  V  ,  MMMMM  ) 

THIS    IS    TABLE    ENTRY    NUMBER    49    OF     LAPLACE    TRANSFORMS     8Y    SPIEGEL 

IMPLICIT    ?E1L*8     (  A, 3,  C-hO-ZJ 

I KPLICIT    C0"PLEX*16     { CJ 

MMMMM     =     M'jIMMM     ¥      1 

ZERO    =    0.0+0 

C  SI    =    DCMPLX(X.Y) 

XPT5    =     X*T 

XPTI     =    y*T 

CXPT    =     CCMPLX(XPTR,XPTI) 

CS3   =    COEX^trx^T) 

CNUM    =    CS1*CS3 

CS4   ■    CS1*CS1*CS1 

CS5    =    0C^PLX(  AL,  ZERO) 

CS6    =    CS5*CS5*CS5 

CDENOM    =    C  S4    +    CS6 

CRESL""    =    CNJM/COENO" 

U    =    DREAL(CRE5L~) 

V  =    DIMAG(CRESLT) 

RF7JRN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.3 

yF  ■  2.2 

MMMMM  =23 

LOSC   =  k 

SUM  =  7.6682  5D+0 

SUMA  =  7.92326D+0 

TIME  =  0.52U16D-1 

ANALYTICAL   =   7.66825D+0 
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Case   50 


£ 


-1 


**  +  a* 


■r     «"•'  +   2e"'i  cos 


•/3aty 


SUBROUTINE    VftLUc    [X tYf T»ALf U,V»MMMMM) 

THIS     IS    T&3LE    ENTRY    NUV6£R    53    OF    LAPLACE    TRANSFORMS    3Y    SPIEGEL 

IMOLI'IT  REAL*8    ( 4  ,3  ,0-H  ,  G-Z) 

IMPLICIT    Ci"PLEX*lb    (C) 

MMMMM    =    mmmmh    +    J. 

ZERO    =     0.3  +  0 

CS1    =    DCMDLXU.Y) 

CS2    =    CS1*CS1 

XPTR    =     x*T 

XPTI     =    Y*T 

CXPT    =     OCMPLXUPTR.XPTI) 

CS3    =     CDEXP(CX^T) 

CN'JM    =    CS2*CS3 

CS4   =    C  S1*C  S1*CS1 

CS5    =    0C*DLX(4L,ZER3) 

CS6    =    CS5*C^5*CS5 

CDENOM    =   "S4   +    CS6 

CRESLT     =    CNL)*VC3EN0"1 

U    =    3REAL(C*ESLT) 

V   =    OIMAG(CRESLT) 

RETJRN 

END 


AL   =   1.25D-1 
T   =   8.D+0 
AA   =   O.k 
Xp  =   -4.89 
yF=  2.3 

MMMMM  =   2U 

LOSC    =   5 

SUM   =   8.34720D-1 

SUMA   =   3.88818D+0 

TIME   =   0.55796D-1 

ANALYTICAL   ■   8.3^720D-1 
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Case   5 


£ 


-I 


«a-o» 


e-"'1  I 

3a*    \ 


rKii/a  _  cos 


v7- 


-    \/3! 


V3 


'I 


SUBROUTINE    VALUE     ( X, Y , T, 4L, U, V , MMMVM) 

THIS     IS    'ABLE    ENT  RY    NUMBER     51    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLi: IT    REAL*8     ( A,B,0-H,O-Z) 

ILLICIT    rnMoL£x*16    I  C  ) 

MfMMM.    =    MMMMM    +•    1 

ZERO    =    O.D+O 

CS1    =    DCMPLX(X,YJ 

XPTR    =    X*T 

XPTI     =     Y*T 

CXPT    =     0C*PLX(X3TR,XPTI  ) 

CMJM   =    CDEXP(CXPT) 

CS4   =    CS1*CS1*CS1 

CS5    =    DCHPLX (AL, ZERO) 

CS6    =    CS5*CS5*CS5 

CDEM1M    =    CS4    -    CS6 

CPESL"     =    CMUM/COEN3M 

U    =    3REALICRESLT) 

V    =    OIMAG(C'RESLT) 

RETURN 

END 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.4 

yF  ■  2.2 

MMMMM   =23 

LOSC    =   3 

SUM   =    3.25349D+1 

SUMA   =   3. 2 5^3^0+1 

TIME  =  0.61204D-1 

ANALYTICAL  =  3-253^90+1 
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Case   52 


£ 


-J 


**-oJ 


3a 


>/3, 


V2at 


V5at 


+  «*>«'* 


SU8R0UTIME    VALUE    (X ,Y ,T , AL,U , V , mmmmm) 

TIS    IS    TABLE    ENTRY    NLM3ER    52    OF    LA°LACE    TRANSFORMS    8Y    SPIEGEL 

IMPLICIT    ">EfL*3    (  A, 3,  3-H,0-Z) 

I  VPLICIT    COMPLEX* 16     (C) 

MMMMM    =    MMMMM    +    1 

ZERO    =     0.0+0 

CS1    =    OCMPLX(X,Y) 

XPT?    =    X*T 

XPTI     =    y*T 

CX°T    =     DCMPLX  (X^^RtXPTI  } 

CS3   =    CDEXP(CX=  T) 

CN'JM    =     CS1*CS3 

CS4   =    CS1*CS1*CS1 

CS5    =     OC^DLX(AL,ZERQ) 

CS6    =    CS5*CS5*CS5 

CDENDM    «   CS4   -    C56 

CRESLT    =    CNU^/COE'lOM 

u  =  dreahcre3l7) 
v  =  dimagcreslt) 
return: 

END 


AL    =    1.25D-1 
T   =   8.D+0 
AA   =    0.4 
Xp    =    -4.89 

yF  =  2.3 

MMMMM   =   24 

LOSC    =   4 

SUM  =  8.33492D+0 

SUMA   =   9.S6283D+0 

TIME   =   0.67418D-1 

ANALYTICAL   =   8.33492D+0 
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Case   53 


Z 


-   *(■ 


cat    +    2«-al/l  cos 


VTai' 


SUB3D 
THIS 
IMDLI 
I  "PL  I 

M  MM  M  M 

ZE«0 
CS1    = 
CS2    = 
XPTR 
XPTI 
CXPT 
CS3    = 
CWH 
CS4    = 
CS5   = 
CS6    = 
C0EN3 
CRESL 
U    =    D 
V   =    D 
RETJR 
END 


C 
D 
C 

M 

T 

RE 

IM 

N 


IME 

TAB 
T   RE 

t  cn 

0.0  + 

a  pl 
si*c 
x*t 

Y*T 

DC  MO 

DEXP 

CS2* 

S1*C 

CMPL 

S5*C 

=  cs 

=  C\j 
AL(C 
AG(C 


V*LUE    (X,YfT,AL,U,V,MMMMM) 

LE    ENTRY    NUMBE'    53    OF    LA°LACE    TRANSFORMS    BY    SPIEGEL 

&L*8    ( A,B  iD-HiO-ZJ 

'">LEX*lt>    (C) 

mm    +    l 

0 

X(X,Y) 

SI 


LX(XPTR,XP"' 
(CXPT) 
CS3 
C1*CS1 

X(  AL, ZERO) 
S5*CS5 
4    -    CS6 
UM/COE  Ol 
PE5LTJ 
RESLT) 


I) 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.3 
Xp  =   -4.99 

yF  =  2.3 

MMMMM  =   2U 
LCSC    =   5 

SUM  =   1.16806D+0 
SUMA   =   2.25373D+0 
TIME   =   0.70200D-1 
ANALYTICAL   =    1.16806D+0 
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Case   5^- 


£ 


-1 


i*  +  4o* 


— 7 (sin  at  cosh  at  —  cosaisinhat) 
4os 


SU6=")UT 
THIS    IS 

i  yp li  :  i 

IVPLIC1 

MMMMM    = 
ZERD    = 
CS1    =    D 
XPTR    = 
XPTI     = 
CXPT    = 
CNUM    s 
CS4   =    C 
CS5    =    0 
CS6    =    4 

C0E\m 

C  RE  S  LT 
U  =    DRE 
V    =    DIM 
RETURN 
END 


IME    VALUE     (X,Y,T,ALtUtV,MMMMM) 

~A8LE    ETRY     N!J*8E3    54    OF    L\0L4CE    TRANSFTRMS    BY    SPIEGEL 
T    RE*L*8    ( A,R,0-H,D-Z) 
T    rnMDi_EX*lo    (C) 

MMMMM    +     1 
O.D+0 

CMOLX  (X,  Y) 
X*T 
Y*T 

coexp(cxpt; 

S1*CS1*CS1*CS1 

CMPLX(  AL.ZERO) 

.D+0*(CS5*CS5*CS5*CS5) 

=    CS4    +   CS  6 

=  CNUM/COENOM 

AL(CPESLT) 

AG(C?ESLT) 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =   0.1* 
Xp   =   -k.kk 

y„  =  2.2 

MMMMM   =   23 

LOSC    =   2 

SUM  =  8.U9272D+1 

SUMA   =   8.U9272D+1 

TIME   =   0.67782D-1 

ANALYTICAL  =    8.49272D+1 
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Case   55. 


X 


-I 


t*  +  4a* 


sin  at  sinh  at 
2a» 


SUBR1U 
THIS     I 
I  MP  LI  : 
IM3L  IC 
MfMVM 
ZER3    = 
CS1    = 
XPTR    = 
XPTI    = 
CXPT    = 
CS3    = 
CNIW    = 
CS4    = 
CS5   = 
CS6    = 
COENOM 
C  RE  SL  T 
U    =    DR 
V    =    DI 
RETURN 
END 


TINE  V4 
S  TABLE 
IT    REAL 

IT    r1^D 

=       KMMMM 

O.D  +  O 
OCMOLX( 
X*~ 
Y*T 

DCMPLX 
COEX  =  (C 

CS1*CS 
CS1*CS1 
DCMPLX( 

*   CS4 

=  CNUM 
EAL(CRE 
MAG(CRE 


LUE     ( Xt Y,T,AL,U,V,MMMMM) 

ENTRY    NIIMBE'     55    IF    LAPLACE    T 
a  o      /  a     a     n_u     ->_  7 \ 


*3     ( At3,0-H,3-ZJ 
LEX*lo    (  C  ) 
+    1 

X,  Y) 


IXPTR, XPTI ) 

XPTJ 

3 

*CS1*CS1 

AL,ZER0) 

CS5*CS5*CS5*CS5) 

+    CSd 

/COENOM 

SLT) 

SLT) 


RANSF3RMS    BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    O.U 

XF    =    -4.5^ 

yF  -  2.2 

MMMMM   =   23 

LOSC    =    3 

SUM  =  3.16447D+1 

SUMA  =  3.16532D+1 
TIME  =  0.5935?D-1 
ANALYTICAL  =  3.l6*V*7D+l 
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Case   56. 


X 


4a* 


■=—  (sin  at  cosh  at  +  cos  at  sinh  at) 


SU8R0U 
THIS  I 
IMdlic 
IPPLIC 

ZER3    = 
CS1    = 
CS2    = 
XPTR    = 
XPTI     = 
CXPT    = 
CS3    = 
CNUM    = 
CS4    = 
CS5    = 
CS6    = 
C DENIM 
CRESLT 
U   =    DP 
V    =    01 
RETURN 
END 


TINE    VALUE     (X,Y,T,ALtU,V,MMMMM) 

S    T*PLE    ENTRY    NU^3E"    56    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IT    ?  EAL*8    (  A,8,0-H,0-Z) 

IT    :0MPLEX*16     ( C) 

=     MVM^M    >     1 

0.0+0 
DCMPLXU  fY) 
CS1*CS1 

X*' 

Y*T 

OC^PLX(  XPTR,  XPTI  ) 
COEX  P(CXPT) 

CS2*CS3 
CS1*^S1*CS1*"S1 
DCMPLX(ALtZERU 
4.D+0*(CS5*CS5*CS5*CS5J 

=    CS4    +    CS6 

=    CNUM/CDENOM 
EAL(CRESL~) 
M  AG(C^ESLT) 


AL    =    1.25D-1 
T   =   8.D+0 
AA   =    0.4 
XF   =   -4.89 

yF  =  2.3 

MMMMM   =24 

LOSC    =   4 

SUM   =   7.73369D+0 

SUMA   =   9.23364D+0 

TIME   =   0.762C6D-1 

ANALYTICAL   =   7.73369D+0 
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Case  57. 


£ 


-1 


$*  +  4o« 


cos  at  cosh  at 


SUBROUTINE    VALJE    <  X  ,  Y,  T,  ALt  U,  V  ,  MMM^m) 

THIS     IS    -46LE    E'lT^Y    NUMBER    37    OF    LA°LACE    TRANSFORMS    BY    SPIEGEL 

IMPLICI T    ' EAL*3    (  A, 3  ,0-H,0-ZJ 

IMPLICIT    COMPLEX*  16    (C) 

MMMMM     =    MMMMM    +    1 

ZERO    =     0.0+0 

CS1    =    DCM°LX{XtYJ 

CS2    =    CS1*CS1*C51 

XPTR    =    X*T 

XPTI    =    y*^ 

CXPT    =    DCMPLX  (XPTR,  XPTI  ) 

CS3   =    COEXP(CXPT) 

CNUM    =    CS2*CS3 

CS4   =    CS1*CS1*CS1*CS1 

CS5   =    OCMPLX( AL.ZERO) 

CS6    =    *.D+0*(CS5*CS5*CS5*CS5J 

CDENOM    =    CS4-    +    CS6 

CRESL'    =    GNU" /CDENOM 

U    =    DREAL(CRESLTJ 

V   =    DIMAG(CRESLT) 

RETJRN 

END 


AL  =    1.25D-1 

T  =   8.D+0 
AA   =    Q.k 
XF   =    -fc.89 

yF  -  2.3 

MMMMM   =   2k 

LOSC    =   5 

SUM  =   8.33730D-1 

SUMA   =   3.85232D+0 

TIME   =   0.82576D-1 

ANALYTICAL   =   8.33730D-1 
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Case  58, 


z 


-1 


i* -a* 


y-y(sinh  at  —  sin  at) 


simouT 

THIS      IS 
1MPLICI 
IMOLICI 
MWMMM    = 
ZEfO    = 
CS1    =     D 
XPTR    = 
XPTI     = 
CXPT    = 
CNUM    = 
CS4  =    C 
CS5    =    D 
CS6    =    C 
CDE^DM 

CRESL' 
U   =    ORE 
V    =    DIM 
RETURN 
END 


INE 

TA 

T    R 

t  r 

MM 

0.0 
CM° 

X*T 

Y*T 

DC^ 

CDE 

SI 

CMP 

S5* 

=  c 
=  c 

Al( 

AG( 


VALUE     (X,Y,T,/\L,U,V,MMMMM) 
BLE    E'JTRY    MUM6E=     55    2?    LAPLiCE    TRANSFORMS    BY    SPIEGEL 
EAL*3     ( A,5,C-h,0-ZJ 
CiMPLEX*16    (CI 
«MM    +    i 

0 
LX(X,YJ 


PLX(  XPTR,  XPTI  ) 
XP(  CXPr  } 

:si*csi-csi 

LX(AL,  ZE3ri) 

CS5*C55*CS5 

S4   -    :S6 

N'UM/CDENOM 

CRESLT) 

CRESLT) 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =    O.k 
X?   =    -k.V* 

yF   =  2.2 

MMMMM   =23 

LOSC    =   2 

SUM   =   8.5^3^90+1 

SUMA   =   8.5^3^905+1 

TIME   =    0.7UU90D-1 
ANALYTICAL   =   8.5^3U9D+1 
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Case   59 


X 


-1 


=-y (cosh  at  —  cos  at) 


THIS 
I  MPLI 
I.M3LI 

MKVMM 
ZERO 
CS1   = 
X  PTR 
XPTI 
CX°T 
CS3    = 
CNUM 
C54    = 
CS5    = 
CS6    = 
CDE^O 
CRESL 
J    =    D 
V    =    D 
RETUR 
END 


UTINE    VALUE     ( X, Yt T.ALt U.V.MMMMM) 

IS    TABLE    EMT^Y     rJM3£R     59    3F    LAPLACE    TRANSFORMS    8Y    SPIEGEL 
CIT    REAL*8     < A,8t J-h,H-Z) 
C IT    CD^LEX^lto    ( C  ) 

=     MM«VM     +•     I 

=    O.D+0 

ocmolxcx,  y) 

=      X*T 

=     Y*T 

=     DCMDLXtXPTR.XPTI  ) 

COEX  P{CXPTJ 
=    CS1*CS3 

CS1*CS1*CS1*CS1 

DCMPLX(AL,ZEROJ 

:  S5*CS5*CS5*CS5 
M    =    CS4    -    CSo 
T    ■   CNUM/COENOM 
'EALIC'ESLTJ 
IMAG(CRESLTJ 
N 


AL   =    1.25D-1 

T   ■   8.D+0 
AA   =    0.3 

y,    =2.2 

MMMMM   =23 

LOSC    =    3 

SUM   =   3.20389D+1 

SUMA   =   3.20897D+1 

TIME   =   0.7S390D-1 

ANALYTICAL   =    3.20889D+1 
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Case   60 


£ 


-1 


i«  -  a* 


27<sinho(  4-  sin  at) 


SUBROUT 
THIS  IS 
IMPLICI 
I MPLICI 

M'.|WMW      = 

ZERO    = 
CS1    =    D 
CS2   =    C 
X  PTR    = 
XPTI     = 
CXPT    = 
CS3    =    C 
CNUM    = 
CS^    =    C 
CS5    =    D 
CS6   =    C 
COENOM 
CRESLT 
U  =    0°E 
V    =    DIM 
RETURN 
END 


INE     VALUE     (X  ,Y  ,T,&L,'J,V  ,W«) 

TA9LE    ENTRY    'JUfBER    60    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 
T    RE4L*8     (  Af  *,  G-H.O-Z  ) 
T    C0MPLEX*16     (CJ 

MMMMM    ♦■     i 
O.D+0 
CMPLX(X,  Y) 
S1*CS1 
X*t 
Y*T 
CCMPLX(X=>TR,XPTI  ) 

DEXO(CXPTJ 
CS2*CS3 

Sl*rSl*CSl*CS  1 
CMPLXULfZEROJ 
S5*CS5*CS5*CS5 
=    CS4    -    C5  6 
=    CNUM/CDENOM 
AL(CRESLT) 
AG(CRESLT) 


AL   =    1.25D-1 
T   =   8.D+0 
AA   ■    0.3 

x     =  -4.54 
y     =  2.2 

MMMMM   =   23 

LOSC    =   4 

SUM  =  8.06669D+0 

SUMA  =  8.32767D+0 

time  =  0.87672D-1 

ANALYTICAL  =  8.06669D-t-0 
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Case  61 


£ 


-1 


$*-  a* 


|  (cosh  a  t  +  cos  at) 


SUBOUT 
THIS     IS 

i  w>  li  :  i 

IMPLICI 

MPMMM    = 

ZERO    = 
CS1    =     D 
CS2    =    C 
XPTR    = 
XPTI    = 
CXPT    = 
CS3   =    C 
CNJM    = 
CS4   =    C 
CS5    =    3 
CS6    =    C 
CDENDM 
CRESLT 
U    =    0"E 
V    =    OIM 
RETURV 
END 


INE     VALUE     (X, Y,T,4L,U, V,MMMMM) 

TABLE    ENrRr    NUMBER    61    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 
T    R5AL*8     ( A,6,D-H,0-Z) 
T    r-IMPLEX*l6    (  C  ) 

yviMMM    +    1 
O.D+O 

CMnLX{X.YJ 
S1*CS1*CS1 
X*T 

DCM°LX(XPTR,XPTI) 

DEXP(CXPT) 

CS2*CS3 

S1*CS1*CS1*CS1 

C^LXt AL .ZERO) 

S5*CS5*CS5*CS5 

=    CS4   -    CS6 

=    CNUM/COENOM 

AL  (CRESLT) 

AG(  CRESLT) 


AL   =    1.25D-1 
T    =    8.D+0 
AA   =    0.3 
XF    =    -4.99 

yF  =  2.3 

MMMMM   =   24 

LOSC   =   5 

SUM  =    1.04169D+0 

SUMA  =  2.13326D+0 

TIME  =  0.9331^0-1 

ANALYTICAL   =    1.04169D+0 
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Case   62 


£ 


-1 


,-bt  _  .-at 


2(6  -  a)y/70 


SUBROU 
THIS  I 
I  MP  LI' 
IMOLIC 
D I  MEMS 
MMMMM 

ZERO    = 
CS1     = 
CS2    = 
CS3    = 
XPTR    = 
XPTI    = 
CXOT    = 
CNUM    = 
CS4   = 
CSS    = 
CALL    C 
CALL    r 
CDENOM 
CRESLT 
U    =    Df 
V    =    DI 
RETURN 
END 


TINE 
S  "A 
IT  R 
IT 
ICN 
=  MM 
0.3 
DC  MP 
DC  MP 
DCMP 

X*T 
Y*t 

DCM 

CDS 

CS1 

CS1 

POWE 

=   C 
MAG< 


VALUE 
PLE  EN 
EAL*8 
O^PLfrX 
MC ALL ( 
MMM  + 
+0 

LX(X,Y 
LX< AL, 
LX (BE, 


(X,     Y,     T,    AL,     BE,    U,    V,     MMMMM,     NCALL) 
T»Y    MUM9E"    62    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 
( A,S,0-H,3-Z) 
*16    (  C  ) 
5  ) 
1 

) 

Z5F0) 
ZERO) 


°LX  (X°TR,XPTIJ 
XP(CXPT) 
♦    CS2 
C33 


+ 

R     (CS4,CS6,5.D-1 ,1 ,NCALLJ 

»     (CS5,CS7,5.D-1,2,MCALL) 

S6    +    C3  7 

NUM /CDENOM 

C<?ESLT) 

CRESLTJ 


AL  =    1.25D-1 
BE   =    -1.25D-1 
T   =   8.D+0 
AA   =    0.3 
Xp    =   -4.52 

yF  =  2.2 

MMBJMM   =   23 

LOSC   =   5 

SUM  =   1.17209D-1 
SUMA   =   6.50880D-1 
TIME   =   0.13301D+0 
ANALYTICAL   =    1.17209D-1 
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Case  63 


X 


-1 


iV«  +  <* 


erf  Vat 


SU5<n 

THIS 

IM°LI 

IfPLl 

DIMEN 

ZERO 

CS1    = 

CS2    = 

XPTR 

XPTI 

CXPT 

CMJM 

CS3   = 

CALL 

CDENO 

CRE3L 

U    =    D 

V    =    D 

RET'JB 

END 


UT!E    VALJE    (Xf     Y,    T,    4L,     U,    V,    MMMMM,    NCALL) 
IS    TABLE    ENTRY     N'JYPER    63    OF    LAPLACE    TPANSF^RMS 
CIT   fl  E-5L-8    (  a,8,D-H,G-Z) 
CIt   Cn^PLEX*16    (  CJ 
SI  ON    NCALU5J 

=    O.D+O 

DCMOLX(  X,Y) 

DC«PLX( AL.ZERD) 
=    X*T 
=     Y*T 

=     OC,JIPLX(XO'R,XPTI  J 
=    CDEXP(CXPT) 

CS1    +    CS2 
C?3'-JEP     (CS3,  CS4,  5.D-1,  l.NCALU 
M    =    CS1*CS4 
T    =    CMUW/CDENOM 
REAL(CRESLT) 
IMAG(CRESLT) 
M 


BY    SPIEGEL 


AL  =    1.25D-1 

T   =   8.D+0 
AA   =    0.2 

yF  =  2-2 

MMMMM  =23 

LOSC   =    5 

SUM   =   2.38352D+0 

SUMA   =   2.64736D+0 

TIME   =    0.85202D-1 

ANALYTICAL   =   2.38352D+0 
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Case  Gk 


£ 


-1 


i 


—         e"  erf  Vat 


SUBROUTINE     VALUE     (X,     Y»    7,    AL»    U,    V,    MM"MMt     NC&LL) 

THIS    IS    TiBLt    ENTRY    NUMBER    64    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    <?EAL*3     ( A, 3, C-h,0-Z) 

IMPLICIT   COMPLEX*  16    (C  J 

DIMENS  ICN    VC  ALL  (5) 

MMHMM    =     M-4MMM    +    i 

ZERO    =     O.O+O 

CSl    =    DCMPLX(X,Y) 

CS2    =    DCMPLX(ALiZtRO) 


XD-R    =     X*T 

XPTI     x     Y*T 

CXPT    =    DC  Mo 

LX(XPTR,XPT 

CNUM    =     COEX 

5( CXPT  J 

CALL    CPCWER 

tCSlfCS3v 

5.0- 

-lf  lvMCALL) 

C  S4    =    CSl    - 

CS2 

COENOM    =    CS3*CS4 

CRESLT    =   CNUM/COENOM 

U    =    OREAL(C 

>.ESLT1 

V   =    DIM  AG(C( 

>ESLT) 

RETURN 

END 

AL  =   1.25D-1 

T  ■   8.D+0 

AA  =   0.4 

X      =    -4.W* 
F 

yF  =  2.2 

MMHMM   =23 

LOSC    =    5 

SUM  =  6.U7907D+O 
SUMA   ■   8.56802D+0 
TIME   =   0.97552D-1 
ANALYTICAL   =   6.47907D+0 
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Case  65 


^^[j^br]  =   "{is-"*  -*»*>} 


SUBROUTINE     VALUE     (X,     Y,    T,     AL ,     BEi    U,     V,     MMMMM,     NCALL) 

THIS     IS    TABLE    ENTRY    NUMBER    65    IF    LAPLACE    TRANSFORMS    8Y    SPIEGEL 


IMPLICIT    REAL*8     (  A,  8,  C-hf  Q-ZJ 

IMOLICIT   C'"JIOLEX*16    (C  ) 

DIMENSION    NCALL  (5) 

MMMMM    =     MMMMM    >    L 

ZE"3    =     0.3+0 

CS1    =    DCMDLX(XfYJ 

CS2    =    OCMDLXt  ALiZEROI 

CS3   =    DCMPLXf  SEtZEROJ 

Xf^R    =    X*t 

XPTI    =    Y*T 

CXPT    =     OrMDLX(X=TR ,XOTI  ) 

CMJM   =    CDEXP(CXPT) 

CS4  •   CS1  -   CS2 

CALL    CPOWEs    (CS4.CS5,  5.0~1»  1, NCALL) 

CDENOM    =   CS3    +    C3  5 

CRESLT    =    C MUM /CDENOM 

U    *   OREAL(CRESLTJ 

V    =    DI  MAG(CRESLT) 

RETJ'M 

END 


AL   =    1.25D-1 

BE  =  3-535533906D-1 

T  =  8.D+0 

AA  ■  0.3 

yF  =  2.2 

MMMMM  =  23 

LOSC  =  5 

SUM  =  1.31286D-1 

SUMA  =  9.058^9D-1 

TIME  =  9.0H6D-2 

ANALYTICAL   =    1.31286D-1 
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Case   66 


£-x 


y/aHT 


J0(al) 


SUBRDU 
THIS     I 

imdli: 

IMPLIC 

DIMENS 

ZERO    = 
CS1    = 
CS2    = 
XP"R    = 
XPTI     = 
CXPT    = 
CNUM   = 
CS3    = 
CALL    C 
CRESLT 
U  =    DP 
V    =    01 
RETURN 
END 


TINE    VALUE     (X,     Y,     T,    *L,     U,    V,    MMMMM,    NCALL) 
S    TABLE    EJT^Y    NUMBER    66    OF    LAPLACE    TRANSFORMS 

IT    REAL*8     <  A,  B,D-H,0-Z) 
IT    C1MPLEX*lo    ( C  ) 
ION    NCALL  (5  ) 
=    ^mmmm    ♦     l 

O.D+0 
DCVPLX(X,Y) 
DCMPLXULtZEROJ 

Y*" 

DCMPLX(  XPTR,  XPTI  ) 

CDEXP(CXPT  ) 
(CS1*CSU    +    (CS2*CS2) 
POWEO     (CSJ,CDESJQM,5.D-1,1  ,  NCALL) 

=   CNUM/CDEMOM 
EAL(CRESLT) 
MAG(CRESLT) 


BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.2 

Xy        =       -U.6U 

yF  -  2-2 

MMMMM   =23 

LOSC    =    5 

SUM  =  7.65198D-1 

SUMA  =  1.22  916D+0 

TIME   =   0.912S6D-1 

ANALYTICAL   =    7.65198D-I 
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Case   67 


£ 


-1 


V«;  - 


/„(»0 


SUSRO 

THIS 

IMPU 

I  MPLI 

DIVEM 

MKMMM 

2E«0 

CS1    = 

CS2    = 

XPTR 

XPTI 

CXPT 

CNUM 

CS3    = 

CALL 

CRESL 

U    =    D 

V   =    0 

RETUR 

END 


UTINE 
IS  TA 
CIT    R 

cue 

S  I  ON 
=    yy 

=     0.0 
DCMP 

DCMP 
=  X*T 
=  Y*T 
=  OCM 
=    CDE 

(CS1 
CPCWE 
1    =   r 

DEAL  < 
I  MAGI 
N 


VALUE     (X,     Y,    T,     4L,     U 

8LE    E'-ITRY    NUMBER    67    3F 
EAL*a    ( A, 3f C-H,0-Z) 

2MDL£X*16      (C) 

MCALL1 5) 

MMV    +•    L 
+  0 

LX(X,  Y) 
LXIAL.ZER3J 


°LX  (X»T3  ,XPT  IJ 

XP(CXPT) 

*CS1)    -    (CS2*CS2) 

R    (CS3,CCENGM,5.D-1,1,MCALL) 

NUM/CDENOM 

CRESLT) 

CRESLTJ 


V,  MMMMMf  NCiLL) 

LAPLACE  TRANSFORMS  BY  SPIEGEL 


AL  =  1.25D-1 
T  =  8.D+0 
AA  =  0.3 
Xp  =  -4.99 

yF  =  2.3 

MMMMM  =  2U 

LOSC  =  5 

SUM  =  1.26607D+0 

SUMA  =  2.31076D+0 

TIME  =  0.10291D+0 

ANALYTICAL   =    1.26607D+0 
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Case  68 


X 


-1 


(y/s'  +  a*  -  a)" 
Vs-  +  a2 


n  >  -1 


a"  ./.(a/) 


SUBROUTINE    VALJE     (X,     Y,     TT     M,     NPOWER,    U,    V,    MMMMM,     NCALL) 
THIS    IS    TABLE    ENTRY    NUMRER    68    np     LAPLACE    TRANSFORMS     BY    SPIEGEL 
IMPLICIT    ?E4L*8    {  A,3f  D-HO-Z) 
I  WPLICIT    CnMPLEX*16    (  C) 
DIMENSION    NCALL  (5) 

MMMMM     =     MM -4  MM     +     1 

ZERO    =    0.0+0 

C  SI    =    OCMPLXf  X,Y) 

CS2    =     DCMP  LX<  AL.ZERO) 

XPT*     =    X*T 

XPTI     =     Y*T 

CXPT    =     CCM'LXUPTR,  XPTI  ) 

CS3    =    (CS1*CS1)    ♦    (CS2-CS2) 

CALL    CPOKER    <CS3,CDENCM, 5.0-1  ,1  , NCALL) 

CS4    =     (CDENQM    -    CSl)**NPnwEP 

CNUM    =    CS4*C0EXO(  CX^T) 

CRESLT    =    CNUM/C0tN0M 

U    =    DREAL<CE$LT) 

V    =    01  MAG(CRESLT) 

RETJRN 

END 


AL  =  1.252-1 
NPOWER  =  3 
T  =  8.D+0 
AA  =  0.4 
Xp  =  -2.16 

yF  =  1.6 

MMMMM  =17 

LOSC  =  1 

SUM   =   3.S2096D-5 

SUMA   =   3.82096D-5 

TIME  =  0.8122UD-1 

ANALYTICAL   =   3-82090D-1 
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Case  69. 


£ 


-1 


(s  -   yV-g2)" 
yfg-  —  a- 


n  >  -1 


J 


a' I, {at) 


SUBPOUTINE    VALUE     (X,     Y,     T,     AL  ,     NDOWER,     U,    V,    Mmmvm,     NCALL) 
THIS    IS    TA6LE    ENTRY    NUMBER    f.9    OF     LAPLACE    TRANSFORMS     BY    SPIEGEL 
IMPLICIT    REU*8    Mf3»0-H,0-Z) 
IMPLICIT    CV'PLEX*lo     (  CJ 
DIME  MSI  ON    NCALL  (5) 
M|MMMV|    =     HMMMM     +     1 

ZERO    =    O.D+O 

C  SI    =    OCMPLXI  X,Y) 

CS2    =    DC^PLXIAL, ZER3) 

XPTR    =    X*' 

XPTI    =    Y*T 

CXPT    =     CC'^PLXU^TR,  X=>TI  ) 

CS3    =    (CS1*CS1)     -    (CS2*CS21 

CALL    CPOWER    (CS3 ,COE NC!M, 5.D-1 , 1 , NCALL) 

CS4    =     (CS1    -    COfcNOMJ**NP"3WE? 

CNUM    =    CS4*CDEXP(CX°T] 

CRESLT    =    :MUM/COEN]M 

U    =    DREAL(CRESLT) 

V    =    DIMAG(CRESLT) 

RETURN 

END 


AL   =    1.25D-1 
N POWER   =   3 
T   =    8.D+0 
AA   =   0.4 
Xp    =   -2.16 

yF  -  i.« 

MMMMM  =17 

LCSC    =    1 

SUM  =   ii.32979D-5 

3UMA   =   ^.32979D-5 

TIME   =   0.7^6720-1 

ANALYTICAL   =   ^.329770-5 
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Case   70 


Z'1 


Vs2  + 


J0(aVt(t  +  2b)) 


SUB ROM 
THIS  I 
IMPLir 

I^PLIC 

01  MENS 

ZERO    = 
CS1    = 
CS2    = 
CS3    = 
XPTt    = 
XPTI    = 
CXPT    = 
CS4   = 
CALL    C 
CS5    = 
CNUM    = 
CRESLT 
U    ■    0* 
V    =    DI 
RETURN 
END 


TINE    VALUE     (X,     Y,    T,    AL  ,     <3E»    U,     V,    mmv^m,    NC4LL) 

S    TA3LE    ENTRY    NUMBER    70    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IT   5  EAL*8    (  i,B  ,0-H,3-Z  ) 

IT    Cr,MPLEX*i6     (C) 

ION    NCALL15J 

=     MMVMM    +     1 

0.0+0 
OCMPLXi  X,  Y) 
DCMPLXJ  M  ,ZEkO  ) 
DCMPLX(S£,ZERO) 

X*t 

Y*T 

DCM0LX(XP7R,XPTI) 
<CS1*CS1)     *     <CSZ*CS2) 
POWER     (CS  *f  COEMOM,  5.0-1,  l.NCALL) 
(CS1    -    CDE>J0,'<)*CS3 

CDEX°(CS5)*C0EX?(CX=>T) 

=    CNUM/CDENOM 
EAL(CRESLT) 
M AG(C^ESLTJ 


AL   =    1.25D-1 
BE   =    -1.25D-1 
T   =   8.D+0 
AA   =    0.2 


Xp    =   -4.64 

yF  =  2.2 


MMMMM   =23 

LOSC   =  5 

SUM   =   7.72088D-1 

SUMA  =   1.23485D+0 

TINE  =   0.3291^0-1 

ANALYTICAL   =    7.72088D-1 
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Case  71 • 


[j0(aVt--b2) 
1° 


t  >  b 
t<  b 


SUBROUTINE    VALUE     (X,     Y,    T,     H,     ^E,    U,    V, 

THIS    IS    TABLE    E'JTRY    NUMBER    71     OF    LAPLACE 

IMPLICIT    REAL*8    (Af3,0-ri,0-Z) 

IfPLICI'    C3MPLEX*Lo     (CJ 

DIMENSION    NCALLI5) 

MMMMM    =     MM vmm    +     1 

ZERO    =    0.0+0 

CS1    =    0CMPLX(  X,  Y) 

CS2    =    DCmdlx  (  AL,  ZERt) 

CS3    =    DCMPLX(-3E,  ZERO) 

XPT?    =     X*T 

XPTI     =     Y*T 

CXPT    =    OCMOLXIXPTR.XPTI  ) 

CS4   =    (CSl*CSi)    ♦    (CS2*CS2) 

CALL    CPOWER    (CS4,CDtNC^,5.D-l, l.MCALL) 

CNUM    ■    CDEXP(C3ENOM*CS3)*CDEXP(CXPT) 

CRESL'     =    CNUM/CDENOM 

U    =   OREAL(CRESLT) 

V    =    DIMAG(C  =  ESLT) 

RETJRN 

ENO 


MMMMM,     NCALL) 
TRANSFORMS    BY    SPIEGEL 


AL  =    1.25D-1 
BE   =    -1.25D-1 
T  =   8.D+0 
AA   =    0.2 

xF  =  -k.6U 
yF  =  2.2 

MMMMM   =   23 

LOSC    =   5 

SUM   =   7.65252D-1 

SUMA   =    1.2U690D+0 

TIME   =   O.86658D-I 

ANALYTICAL   =   7.65252D-1 
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Case   72 


£ 


-1 


i 

(sl+a2)U2 


—  tJi(at) 


SUBROUT 

THIS    IS 
IMPLICI 

I"PLICI 
DIMEMSI 
MMMMy     = 
ZERO     = 
XPTR    = 
XPXI     = 
CXPT    = 
CNUM    = 
XREAL    = 
XIMAG    = 
CS1    =     D 
CS2    =    C 
CALL    CD 
CRESLT 
U   =    ORE 
V    =    JIM 
RETURM 
END 


IME    VALUE    (X,     Y,    T,     \\_,     (j ,    v,    MMMMM,    NCALL) 

T4BLE    ENTRY    NUMBER    72    OF    LAPLACE   TRANSFORMS     BY    SPIEGEL 
T    R  E1L*8    {  A, 5  ,D-h,  O-Z  ) 
T    CT*DLEX*i6     (  C) 
ON    NCALU  5) 

o.o+o 

X*T 
Y*T 

OCMPLX  (XPTR,XPTI) 
COE  XP(CXPT) 
(X*X)     -    (Y*Y  )     +    <AL*4L  ) 

2.o  +  o*x*y 
gmplx(x*eal»ximag) 

S1*C<1*CS1 

OHER    (CS2,    COENCP,    5.3-1,     1,     NCALL) 

=    CMU'VCOENOM 
AL(CRESLT) 
AG<  CRESLT) 


AL   =    1.2  5D-1 
T  =   8.D+0 
AA   =   0.4 
Xp    =    -4.44 

yF  =2.2 

MMMMM   =23 

LOSC    =   3 

SUM   =   2.81632D+1 

SUMA    =   2.Q1715D+1 

TIME   =   0.11476D+0 

ANALYTICAL   =    2.81632D+1 
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Case  73. 


o6 


-i 


(,2  +  „2)3/2 


tJ0(at) 


SUBROUTINE    VALUE     (X,    Y,    T,    AL.    U,    V,    MMM»<Mf    MC4LL) 

ThIS     IS    T4PLE    ENTRY    NUMBER    73    OF     LAPLACE    TRANSFORM 

IMPLICIT    REAL*3     (  A  ,  8  ,  C-  F,  O-Z  ) 

IMPLICIT  COMPLEX*  16    (  C) 

DI^ENS  ION    NCALLt  5) 

MPVMM    =     MMWMM    *■    1 

ZER3    =     C.D-t-0 

CS1    =    DCMPLXiXtYJ 

XPTR     =     X*T 

XPTI     =     Y*T 

CXPT    =     DCMPLX  (X5TP,XPTI) 

CS2    =    CDEXP(CXPT) 

CNUM    =     CS1*CS2 

XREAL    =     (X*X)     -     (Y*Y)     +    (AL*fiL) 

XIM*G    =    2.3+0*X*Y 

CS3    =     DCMOLX(XRElL,XIMAG) 

CS4    =    CS3*CS3*CS3 

C4LL    CPOUER     (CS4f    C3ENCV,    5.3-1  f     1,    NCALL) 

CRESLT    =    CNU^/CDENOM 

U   -    DREAL(CRESLT) 

V    =    3IMAG(CRESLT) 

RETJRN 

END 


S    BY    SPIEGEL 


AL  =    1.25D-1 

T   =   8.D+0 
AA   ■    0.4 
XF   =   -4.89 

YF   =   2.3 

MMMMM   =   ZUr 

LOSC    =   4 

SUM  =  6.12158D+0 

SUMA    =   8.3^9983+0 

TIME   =    0. 12620D+0 

ANALYTICAL   =   6.12158D+0 
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Case  7k.  t 

y'1 


(•*  +  a-')3"  


SUBROUTINE    VALUE     (X,     Y,    t,     «l,     u,     V,    MMMMM,    NCALL) 

THIS    IS    TABLE    E  1TRY    NUMBER    74    OF     LAPLACE    TPANSFOR^S    BY    SPIEGEL 

IMPLICIT    RE*L*3     (  A,  3  ,  D-h,3-Z  ) 

I«*°LIC  IT    CT*°LEX*16     (CJ 

DIME^S  ION    MC»LL< 5  J 

MMMMM    =     M^MMM    +    1 

ZER3    =    O.O+O 

CS1    =    DCMDLXU.Y) 

XPTR     =    X*- 

XPTI     =     Y*T 

CXPT    =     DO*PLX<XPTR, XPTI ) 

CS2    =    COEXP(CXPT) 

CS3  =  csi*rsi 

CNUM    =     CS2*CS3 

XREAL    =    (X*X)    -    ( Y*v)     +     (AL*AL) 

XIMAG    =    2.D+0*X*Y 

CS4    =    DC^PLX(XREAL,X IMAGJ 

C  S5    =    C  S4*CS4*C  S4 

CALL    CPOWER     (CS5,    CDENOM,     5.3-1,      1,     NCALL ) 

CRESL"    =   r.NUM/CDENOM 

U   =    3REALICRESLTJ 

V    =    DIM  AG(CR<ESLT) 

RETURN 

END 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =0.4 
Xp    =   -4.89 

y¥  =  2-3 

MMMMM  =   2U 

LOSC    «    5 

SUM   =    3.251U7D-1 
SUMA   =    3.72  5^70+0 
TIME   =    0.107720+0 
ANALYTICAL   =    3.251>7D-1 
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Case   75. 


£ 


-1 


(»=  -  a=)3 


tl,(at) 


SU350U 
THIS     I 

impli: 

I^PLIC 
0IMEN5 
MMM.MM 
ZERO    = 
XPT*     = 
XPTI     = 
CX°t    = 
CNUM    = 
XPEiL 
XI  MAG 
CS1    = 
CS2   = 
CALL    C 
CPESLT 
•J    =    0° 
V   =    DI 
RETURN 
END 


TIME    VALUE     (X,      Y,     T,    4L,     Ut    V,    MMMMM,     NCALL) 

S    X4BLE    E4TRY    DUMBER    75    ">F    LAPLACE    TRANSFORMS    BY    SPIEGEL 

IT    REAL*8    ( A,3,0-H,0-Z) 

IT    COMPLEX* 16    IC  ) 

ION    NCALL  (5  ) 
=    M"MMM    •►     1 

0.0+0 

X** 

DCMPLXCXPTR,  X°TI  ) 

COE  XP(CXPT) 
=     <X*X)    -    (  Y*Y)    -     (AL*4L) 
=    2.0+0*X*Y 
DCMPL  X(  XREAL.XIMAG) 
CS1*CS1*CS1 
POWER    (CS2,     COENOM,    5.0-1,     1,     NCALL) 

=   CNUM/CDENCM 
EALICRESLTJ 
MAG(CRESLTJ 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.5 
Xp  =  -4.3^ 

yF  =  2-z 

MKMMM  =23 

LOSC    =    3 

SUM   =   3.617C2D+1 

SUMA   =    3.638U6D+1 

TIME  =  0.11089D+0 

ANALYTICAL   =    3.6l?02D+l 
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Case  76 


V?  (,.'  _  a»)3/4 


tl0(at) 


subrou 

THIS    I 

impli; 

I*PLIC 

D  I  WENS 

ZERO    = 
CS1    = 
XPTR    = 
XPTI    = 
CXPT   = 
CS2    = 
CNUM    = 
XRE4L 
X  I  "AG 
CS3    = 
CS4   = 
CALL    C 
CRESLT 
U   =    OR 
V    =    01 
RETURN 
END 


T  INE 
S  TA 
IT  R 
IT  C 
I  ON 
=     MM 

0.0 
OCMD 

X*T 

DCM 
CDEX 

CS1 
=  (X 
=  2. 
DCMP 
CS3* 
POWE 

e\1 

MAG( 


V£L'JE     (X,     Y,    T,     Mi     U,     V,    MMMMM,    NCALL) 
BLE    EN~RY     NUMBER    76    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 
EAL*8     ( A,3t0-H,0-ZJ 
0MPLEX*16     (C) 
NCALL (5; 
mmm    +    1 

+0 

lx(x,  yj 


DLXIXPTR, XPTI ) 

C(CXCT  ) 

*CS2 

*X)     -     (Y*Y)     -     UL*AL) 

D+0*X*Y 

LX(XREALtXIMAG) 

CS3*CS3 

3     (CS*,     COENOM,    5.0-1, 

NU^/COENOM 

CCESL') 

CRESL") 


1,     NCALL) 


AL   =    1.25D-1 

T  =   8.D+0 
AA   =   0.5 
Xp   =    -^.79 

yF  -  2.3 

MMMMM   =   2b 
L05C    =   b 

SUM   =    1.01285D+1 
5UMA    =    1.52218D+1 
TIME    =    0.13250D+0 
ANALYTICAL  =    1.012S5D+1 
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Case  77 


Y>      '•*■       (7a  J~a:)—>  ;o(«»0   +  atl{(at) 


SUBROUTINE    VALUE     (X,     Y,    T,     AL,    U,    V,     MMMMM,     NCALL) 

THIS    IS    TABLE    ENTRY    NUMBER    77    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT    3EAL*8     (  A, 3,  0-HtO-Z) 

IMPLICIT    C0'1PLEX*16     (CJ 

DIMENSION    MT.  ALL  (5) 

MMMMM     =     MMMMM    +    1 

ZERO    =    O.D+O 

CSl   =    OCMDLXt  X,  Y) 

X  P'R    =     X*T 

XPTI     =    Y*T 

CXPT    =     OC"DLX(XPTP  ,XPTI  ) 

CS2    =    CDEXPICAP1") 

C  S3   =    CS1*CS1 

CNUM    =    CS2*rS3 

XREAL    =     (X*X)    -    <Y*Y)     -     (AL*AL) 

XIM4G    ■    2.D+0*X~Y 

CS4    =    DCMPLX(XREAL»XIMASJ 

C  S5    =    CS4-*  C  C4*CS4 

CALL    CPOWEP~(CS5,     CDENQM,    5.3-1,     1,    NCALL) 

CRESLT    =    CNUM/CDEN3* 

U   =    DREAL(CDESLT) 

V    =     DIM  AG(C°.5SLT) 

RETURN 

EM 


AL   =    1.25D-1 

T  =   8.D+0 
AA    =    0.5 

Xp    =    -4.79 

yF   =  2-3 

MMMMM  =   24 

LOSC    =   5 

SUM  =    1.831223+0 

SUMA    =    7.71705C+0 

TIME   =    0.11807D+0 

ANALYTICAL   =    1.83122D+0 
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Case  78 


£ 


-1 


SU8R0UT 
THIS  IS 
IMPLICI 
I ^PLIC  I 
OIMENSI 

MMMMM      = 

ZERO    = 

CSl    =    D 

CS2    =    D 

CS3    =    - 

XP-R 

XPTI 

CXPT 

CNUM 

CALL 

CR£SLT 

U    =    ORE 

V    =    DIM 

RETURN 

END 


CP 


I  ME 

TA 
T  P- 
T    C 

ON 

MM 

0.0 

CMC 

CM  0 

CS2 

X*T 
Y*T 
DCM 
CD6 
OWE 
=  C 
AL( 
AG< 


V  ALU 
BLE  E 
EAL*8 
OMPLE 
NC4LL 

MMM  + 
+0 

LX(  X, 
LX(  AL 
/CSl 


E    (X, 
NTRY 

(  A,°. 
X*15 
(5) 

1 


Y,     T-    M 

NUMBER  81 
.Q-H.O-Z) 
(C) 


u 

OF 


V,     MMMMM,     NCALL) 
LA°LACE    TRANSFORMS 


BY    SPIEGEL 


Y) 

.ZERO) 


°LX  (XPTR.XOTI) 

XP(CSj)  *CDEXP(CXPT) 

R    ICS  l.CDENOM, 5.0-1,  1.NCALLJ 

NUM/COENOM 

CCESLT) 

CRESLT) 


AL   =    1.25D-1 
T   =    8.D+0 
AA   =    0.2 

JF 

MMMMM   =23 

LOSC    =    5 

SUM  =  -8.3.0O93D-2 
SUMA  =  ^.^719^D-1 
TIME  =  0.96l^8D-l 
ANALYTICAL   =    -8.30093D-2 


170 


Case  79 


£ 


-1 


gin  2\/o7 

""""  V7a 


SLB*OU 
THIS  I 
I^PLIC 
IMPLIC 
DIMENS 
MMMMM 

ZERO    = 
CS1    = 
CS2   = 
CS3    = 
XPTR    = 
XPTI    = 
CXPT    = 
CNUM   = 
CS4   = 
CALL    C 
CRESLT 
U    =    DO 
V    =    01 
R  E  TUR  N 
END 


TINE 
S  Tt 
I  T    R 

IT    C 

ION 

=     MM 

O.D 

DCMP 

DC  Mo 

-CS2 

X*T 

Y*T 

oc« 

CDE 
CS1* 

POWE 

EAL7 
MAG( 


VAL 
3LE 
EAL* 
OMp|_ 
MCAL 
MMM 
+0 

LX(X 
LX(A 
/CS1 


'JE     (X,     Y,    T,    ALf     U,    V,     MMMMM,     NCALL) 

ENTRY    NUMBER     82    OF    LAPLACE    TRANSFORMS    BY    SPIEGEL 

8     (AtS,C-h,a-ZJ 

EX*16    (C) 

L  (5  J 

*    1 

,y) 

L,ZERJ) 


°LX  IXPT  R,X°TI  ) 

XP(CS3)  »COEXP(CXP") 

CS1*CS1 

R    CCS4, CCENGM,5. 0-1,  1,  NCALL) 

NUM/CDENOM 

C°ESLT) 

CRESLX) 


AL   =    1.25D-1 

T   =   8.D+0 

AA   =    0.2 

X^   =  -4.64 
F 

y„  =  2.2 

MMMMM  =23 

LOSC    =    5 

SUM   =    1.45103D+0 

SUMA  =   1.82676D+0 

TIME   =    0.10018D+0 

ANALYTICAL    =    1.4<103D+0 
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Case  80 


Z 


-1 


r-«" 


»>-l 


—        (£)"".U2^) 


SUBPOUTIME    VALUE    (X,    Y,    T,    AL. 
THIS    IS    TABLE    ENTRY    NUMBER    83 
I^PLIC  IT    o  E4L*8     I  *t3  ,0-HtO-Z  J 
I  MPLICIT   C0vp|_EX*16     (C) 
DIMENSION    MCALL  (5) 
MMMMM    =    mmmmm    +    i 

ZERO    =    O.D+0 
CS1   =    DCMPLX(  X,  Y) 
CS2    =    OCMPLX  (AL,  ZERO) 
XPTR    =    X*T 
X»TI    =     y*T 

CXPT    =     0CM3LX(X5TP,XP"I  ) 
CNUM    =    COEXP(  CS3)*CDEXP(CXP') 
COENTM    =    CS1**(^P0WEP    -►    1) 
CRESLT    =    CNUM/CDENOM 
U   =    DPEAL(CRESLT) 
V   =    DIMAGCCRESLT] 
RETURN 
END 


KjorWP"3,     |Jt     Vt     MMMMM, 

OF     LAPLACE    TRANSFORMS 


NCALL) 

BY    SPIEGEL 


AL   =    1.25D-1 

NPOWER    =    3 
T  =   8.D+0 
AA  =   0.3 
XF   =    ~U.11 

yF  =  2.1 

MMMMM   =   22 

LOSC    =   k 

SUM  =   6.60189D+1 

SUMA   =   6.60189D+1 

TIME  =  6.9758D-1 

ANALYTICAL   =    6.60139D+1 
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Case  81 


r'M  = 


e-a'/At 

•J7i 


SU8R0UTIME    VALUE    (X,     Y,    T,     XI,     ij  ,    v,    MMMM.M,    NCALL) 

THIS    IS    TABLE    E  J7RY    NUM9ES    34-    OF    LA°LACE    TRANSFORMS     BY    SPIEGEL 

IMPLICIT    °EAL*8    (  A,8,0-H,3-Z) 

I MPLICIT   C0MPLEX*16     ( C) 

DIMENSION    ^CALL( 5) 

MMMMM    =     T4M.MM    +•     1 

ZERO    =    O.D+0 
CS1    =    DCMt>LX(X,Y) 
CS2    =    OCMPLX(ALf  ZEROJ 
XPTR    =    X*T 
XPTI     =     Y*T 

CXPT    =     OCMPLX  (X°TR,XPTI) 
CALL    CPOWFR    (CS1 ,CDENCM, 5.0-1 ,1 tMCILL) 
:  S3    =    -CS2*CDEVMH 
CNUM    =    C0EXPICS3}*C0EX?(CXPT] 
CRESLT    =   CMUM/COENOM 
If    =    D^EALCC^ESLT) 
V    =    OIMAG(C=!ESLTJ 
RETUPN 
END 


AL   =   1.25D-1 
T   =   8.D+0 
AA   =   0.2 
Xp  =  -4.64 

yF  =  2.2 

MMMMM   =   23 

LOSC   =  5 

SUM  ■   1.99374D-1 

SUMA   =  6.27^4D-1 

TIME   =   0.93^700-1 

ANALYTICAL   =    1.99374D-1 
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Case   82. 


£ 


.v7 


—      2\/7r* 


,-oV4l 


SUBROUT 
THIS  IS 
IMOLICI 
I  fPLICI 
DIMENSI 
MMMM^     = 

ZERO 
CS1    =    0 
CS2    =    0 
XPTR    = 
XPTI    = 
CXPT    = 
CALL    CP 
CS4   =    - 
CPESLT 
U   =    DRE 
V    =    DIM 
RE"JRN 
EM) 


IME 
TA 

T    '. 

T    C 

ON 
MM 

CMP 

CMO 
X*T 
Y*T 

DCM 

owe 

CS2 
=  C 
AL( 
4G( 


VALJE     (X,     Y,     T,     4L,     Uf    V,    MMMMM,    NCALL) 
3LE    ENTRY    NUM3E?    35    ^F    UPLACE    TRANSFORMS    BY    SPIEGEL 
EiL*8    ( A, 8  ,D-H,0-Z) 
T^PLEX*Lb     (Ci 
NCALL(51 

M«M     +     1 
O.D+O 
LX(  X,Y) 
LXUL,Z5RO) 


PLX(XPTR, XPTI ) 

R    (CS1,C3  3,5.C-1,1,NCALLJ 

*CS3 

DEXP(CS^)*COEXP(CXPT  ) 
CRESLT) 
CRESLT) 


AL   =   1.25D-1 
T  =   8.D+0 
AA   =   0.2 

Xp  =  -5.09 

yF  =  2.3 

MMMMM   =   2b 

LOSC   =   6 

SUM  =   1.5576 ID- 3 

SUMA   =   3.03773D-1 

TIME   =   0.10197D+0 

ANALYTICAL   =    1.55761D-3 
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Case  83 


-1 


t    \_ 


-.si 


erf  (a/2 yft) 


SUBROUTINE     VALUE     (X,     Y,    T,    AL,    U,    V,     MMMMM,     NCALL) 
THIS     IS    -5BLE    ENTRY    NUMBER     66    OF     LAPLACE     TRANSFORMS 
IMPLICIT    REAL*S     I AtBtO-HtO-ZJ 
IMPLICI  T    COMPLEX-  16     (C) 
DIMENS ION    NCALL (5) 
MMMMM    =     MMMMM    ♦■     1 
ZERO    =     O.O  +  O 
ONE    =    1.0+0 
CS1    =    DC«-'PLX(X,  Y) 
CS2    =    OC*°LX<  AL  .ZERO) 
CS3    =    DCMPLX(ONE, ZERO) 
XPTR    =    X*T 
XP+I     =     Y*T 

CXPT    =     OCMPLX  (X°TR,XPTI) 
CALL    CPOWER    (CS1.CS4, 5.0-1  fl  , NCALL) 
CNUM    =     (CS3    -    CDEXP(-CS2*CS4J  )  *COE  XP  (CXPT) 
CDENOM    =    CS1 
CRESLT    =   CNUM/C3EN0M 
U    =    DREAL(C'ESLT) 
V   =    DIMAG(CRESLT) 
RETURN 
ENO 


BY    SPIEGEL 


AL   =    1.25D-1 
T  =   8.D+0 
AA   =    0.2 
Xp   =  -k.6k 
yF  =  2.2 

MMMMM   =   23 

LCSC    =    5 

SUM   =   2.49298D-2 

SUMA   =    8.07480D-2 

TIME   =   0.93158D-1 

ANALYTICAL  =   2.<J-9292D-2 
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Case  84. 


z'X- 


-.vT 


erfc  (o/2-/F) 


SUBRDU 
THIS  I 
I"-10LIC 

IMDLIC 

DIMENS 

MMMMM 

CS1  = 
C  S2  = 
XPTR  = 
XPTI  = 
CXPT  = 
CALL  C 
CNUM  = 
CDENOM 
CRESLT 
U  =  DR 
V  =  DI 
RETURN 
ENO 


TINE  W 
S  'A8LE 
IT  R  EAL 
[t    rnwp 

ion'nca 

ZER->    = 

DCMPLX( 

OCMDLX( 

X*" 

Y*T 

OCMPLX 
POWE"     ( 

CDEXP{ 

=  r.si 

=  CMUV 
EAL  (CPE 
M  AG(CRE 


LUE 

ENT 

*a  ( 

LEX* 
LL(5 
+  1 
0.0+ 
X,Y) 
AL,Z 


(X,     Y,    T,    4Lt     U,    V,    MMMMM,     NCALL) 
RY    TJM3EP    87    IF    LAPLACE    TRANSFORMS 
A  ,3  O-rl,0-Z) 
lb    (CJ 
i 

0 

ERO) 


BY    SPIEGEL 


(XPTR, XPTI) 

C51,  CS4,  5.D-1,  L,NCALLJ 

-CS2*CS4)*CDEXP(CXPT) 

/COENQM 

SLT) 

SLT) 


AL   =    1.25D-1 
T   ■   8.D+0 
AA   =    0.2 

xF  =  -4.6U 

yF  ■  2.2 

MMMMM   =   23 

LOSC    =    5 

SUM   ■    9.75070D-1 

SUMA   =    1.35850D+0 

TIME   =   0.99658D-1 

ANALYTICAL  =   9.75C71D-1 


176 


Case  85 


-7  T    *"*  " 


—  V  ZyftJ 


SUBROUTINE    VALUE    (X,     Y,    %    AL  ,     RE,    U,    V,     MMVMM,     NCALL  ) 

THIS    IS    TABLE    ENTRY    NU«3ER    88    OF     LAPLACE    TRANSFORMS    BY    SPIEGEL 

IMPLICIT     RE0L*8     (  A,  3,  C-H,0-Z  i 

I MPLICIT   C0MPLEX*16     (CJ 

DIMENSION    "J  CALL  (  5  J 

MMMMM    =     MMMMM    +    1 

ZERO    =    0.  0+0 
CS1   =    DCM?LX(X,  Y) 
CS2    =    DCMPLX (AL, ZERO) 
XPT\    =     X*T 
XPTI    =     Y*T 

CXPT    =    DCMt>LX(X°TR,X3TI  ) 
CS3   =    DCMPLXJ3E  ,ZERO) 
CALL    CPOWER     ICSltCS*.  5.0-1,  1, NCALL) 
CNUM   =    CDEXPC-C32*CS4)*C0EXP(CXPT) 
CS5   =    CS3    +    CS4 
COENOM    =    CS4*CS5 
CRESLT    =   CNUH/COENOM 
U  =    DREAKCRESLT] 
V    =    DIMAG(CPESLT) 
RETURN 
END 


AL   =    1.25D-1 

BE   »    3.^57^089060-1 

T  =  8.D+0 

AA   =   0.2 

Xp   =   -4.6U 

yF  ■  2.2 

MMMMM   =   23 

LOSC    =   5 

SUM  =  U.27375D-1 

SUMA  =  7.5575D-1 

TIME   =   9-929^D-2 

ANALYTICAL   «  U.27375D-1 
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Case   86 


£ 


-1 


-(Sf) 


«-»•  -  •-•« 


SUB*} 

ThIS 

IMPLI 

IMPLI 

DIMEN 

MrnnM 

ZERD 
CS1    = 
CS2   = 
CS3   = 
XPTR 
XPTI 
CXPT 
CS4    = 
CS5   = 
CS6    = 
CALL 
CRESL 
U    =    0 
V   =    0 
RETUR 
END 


UTI'IE 
IS    TiR 
CIT    RE 
CI  T   CO 
SIQN    N 

=      *MM 

=     0.0  + 

DCMOL 

OCMPL 

DCMPL 

=    X*T 

=     Y*T 

=     DC  MP 

CS1    + 

CS1    ♦ 

CS4/C 

CLOG    ( 

T     -■   CS 

REAL(C 

IMAGIC 

N 


VALUE 
LE  ENT 
AL*8  ( 
MPLEX* 
CALL(5 
MM  ♦■  1 
0 

X(X,  Y) 
X(AL,Z 
X(6E,Z 


(X,Y,T,AL,8E ,U,V,MMMMM,NCALL ,MCALL  ) 

RY    NUMBER    90    ">F     LAPLACE     TRANSFORMS    8Y    SPIEGEL 

A»8tO-H,0-ZJ 

lo    (  C) 

J,     MCALH5) 


ERD) 
ERQ) 


LX<  XPTR,  XPTI  ) 

CS2 

CS3 
S5 

CS6,     CS7,     1,    MCALLJ 
7*CDEXP(CXPTJ 
RESLT) 
RESLT) 


AL   "    1.25D-1 

BE   ■   2.5D-1 

T  ■  8.D+0 

AA   =    0.2 

Xp   =   -*.21 

yF  »  2.1 

MMMMM   =22 

LOSC    =    5 

SUM  *    -2.90680D-2 

SUMA   ■   1.017*00-1 

TIME   =   7.  009*0-2 

ANALYTICAL   =   -2.906fiOD-2 


178 


Case  87 


£ 


-1 


2* 


Ci  (at) 


SUB'TUTIME    VALUE     (X,     Y,    T,    M,     U, 

THIS    IS    ""ABLE    ENTRY    NUM8E"    91    HF 

IMPLI"  I T    PE4L*3     ( A, 3, D-H.I-ZJ 

Ivt>LICIT    C0'«PL£X*16    (C) 

DIMENSION    MCALL     (5) 

MM^MM    =     MMMMM    ••■     1 

ZERO   =     O.D+O 

CS1    =    DCMPLX(X«YJ 

CS2    =    OCMDLX(  AL.ZERO) 

CS3    =    l(C*l*CSl)    +    (CS2*CS2))/CCS2*CS2) 

XPT>    =    X*T 

X°TI    =     Y*T 

CXPT    =     DC^PLX  (XPTR,XPT  IJ 

CALL    SLOG     (CS3,    C  $4 ,     1,    "CALL) 

CNLM    =     CS4*CDEXP<  CXPT  ) 

CS5    =    DC"PLX<2. O+O, ZERO) 

CRESLT    =   CNUH/<CS5*C$1] 

U    =   OREAL  (CRESLTJ 

V   «    OIMAG(C«ESLT) 

RETJ"N 

END 


Vi     MMMMM,    MCALL) 
LAPLACE    TRANSFORMS    3Y    SPIEGEL 


AL  «  1.25D-1 
T  ■  8.D+0 
AA  =  0.2 

xF  =  -5.09 
yF  ■  2.3 

MMMMM  =  2U 

LOSC  =  5 

SUM  =  -3.37^0^0-1 

SUMA  »  1.27603D+0 

TIME  «  8.82U4D-2 

ANALYTICAL   =    -3O7^04D-l 
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Case   88 


X 


-1 


In  [(*  +  a)/a\ 


Ei(at) 


SU8R0UT 
THIS    IS 
IMPLICI 
I^PLICI 
OIMEMSI 
MHMMM    = 
ZERO    = 
CS1    =     D 
CS2    =    D 
CS3    =    < 
XPTR    = 
XPTI     = 
CXPT 
CALL 
CNUM 
C  RESIT 
U    =    DPE 
V    =    DIM 
RETURN 
ENO 


CL 


INE    VALUE     (X,    Y,    T,     AL,     U 

TA3LE    ENTRY    NUMBER    92    OF 
T    B  EAL*8     (  A,  3,  D-H,  1-Z) 
T    COMPLEX* 16     I CJ 
ON    '•'CALL     (  5) 

MMMMM     ♦     1 

c.o+o 

C"OLX(X,  Y) 

CMD|_X(AL, ZERO) 

CSl    +    CS2J  /CS2 

X*T 

Y*T 

OCMPLX(XPTR,XPTI  ) 

OG    (CS3,    CS4,     1,     'CiLD 

CS4*C0tXP(  CX°T) 

=   CNUM/CS1 

ALCCRESLT) 

AG(CRESLT) 


V,     MMMMM,     MCALL) 
LAPLACE    TRANSFORMS     BY    SPIEGEL 


AL  =  1.25D-1 
T  =  8.D+0 
AA  «  0.2 

xF  =  -5.09 
yF  c  2.3 

MMMMM  =  2U 

LOSC  »  5 

SUM  -  2.193^^0-1 

SUMA  *  i.5509^r>o 

TIME   «=   7.212UD-2 

ANALYTICAL    «   2.1938^-1 
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Case  89 


£ 


-1 


In  ' 


-  (In  t  +  y) 
y  =  Euler's  constant  =  .5772156. 


SUBROUTINE    VALUE     (X,     Y,    T,    AL ,    U, 
THIS    IS    TABLE    ENTRY    NUMBER    96    OF 
ILLICIT   5E'1*8    (  A,3,0-H,0-Z) 
IMPLICIT    COMPLEX*16     (CJ 
DIMENSION    "CALL    (5) 

MHMMM     -      MVMMM     ♦      J. 

ZERO    =    0.0+0 

CS1    =    0CMPLX(  X,  Y) 

CALL    CLOG    (CS1,    CS2,     1,     MC4LL) 

XPTR    =    x*~ 

XPTI     =     Y*T 

CXPT    =     OC^DLX(XPTR,XPTI ) 

CNUM    =    CS2*CDtXP(  CXPT  J 

CRESLT    =   CNUM/CS1 

U    =    OREAL(C°ESLT) 

V    =    01  MAG(CRESLTJ 

RETURN 

ENO 


Vf  MMMMM,  MCALL) 
LAPLACE  TRANSFORMS 


BY    SPIEGEL 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.2 

xF  «=  -5.09 

yF  -  2.3 

MMMMM   =   Zh 

LOSC    =   5 

SUM  ■   -2.6566(fD+0 

SUMA   ■   2.80606D+0 

TIME   '   7.7818D-2 

ANALYTICAL   »    -2.65666D+O 
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Case  90. 


X 


-1 


(Inl  +  r)*   -    Jr* 
y  =  Kuli-r's  conaUnt  =  .577215G. 


SUB'OUT 
THIS     IS 
IMPLICI 
IMPL  IC  I 
0IMEN5I 
MMMMM    - 
ZERD    = 
CS1    *    0 
XP"R    = 
XPTI    = 
CXPT 
CALL 
CNUM    = 
C  RE  SL  T 
U    =    DPE 
V    =    DIM 
RETURN 
END 


CL 


INE     VALUE    (X,     Y,    T,    AL,     U 
TABLE    E^try    NUMBER    97    HF 

T    RfEAL*8    ( A,8  ,0-H.O-ZJ 

T    cnMoL£X*16    (C  ) 

ON    MCALL    (5) 
MMMMM    ♦•      1 

O.D+0 

CMPLX(X,Y) 

X*T 

Y*T 

DCMPLX(XPTR,XPTI  J 

^G    (CS1,    CS2,     It    "CALL) 
CS2*C?2*CD£XP(CXdt) 

=  CNUM/CSl 
AL(CPESLT) 
AG(CRESLT) 


V,    MMMMy,    MCALL) 
LAPLACE    TRANSFORMS 


BY    SPIEGEL 


AL   =    1.25D-1 
T  *   8.D+0 
AA   «   0.2 


XF  -   -5.09 

yF  ■  2.3 


MMMMM    =    2U. 

LOSC    «=   U 

SUM   ■    5.^1289D+0 

SUMA   =    5.1*3^70+0 

TIME   =   7.7376D-2 

ANALYTICAL   =    5. ^12890+0 
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Case  91 


X 


-1 


—  -  f  Y    (at) 


SUBPHlFINE    VAL'.JE    (X  ,Y  ,T  t  AL  ,BE  ,  U  ,  V  ,  Mmmmm  ,  vjr  ALL  t  MC  AL  L) 
THIS    IS    TABLE    ENTRY    *21     PAGE    265    EGBERTS    AND    KAUFMAN 
IMPLICIT    =EiL*8    { A,B,D-H,1-Z) 
I^PLICI  -   COMPLEX*  lb     I  C) 
DIMENSION    MCALL(5)i     *CALL(5) 

MMMMM     =      MM^MM     +     J. 

ZERO    =    O.D  +  0 

C  SI    =     0CMDLX(X,  Y) 

CS2    =    DCMDLX(AL,ZER">) 

XPTR    =    x*T 

XPTI    =     Y*T 

CXPT    =     DC^OLX  (X^TR,XPT  I) 

CS3   =    (<CS1*CS1)    *    <CS2*CS2>> 

CALL    CPHWEC     (CS3tCDEN0Mf 5.3-lf 1,NCALL) 

CS4    =     (CDENOM    +    CSD/CS2 

CALL    CLOG    (CS4,    CS5,     1,     "CALL) 

CNUV    =    CS5*CDEX=(  CX°T  J 

CRESLT    =   CNUM/COENO^ 

U   =    ORE  AL(  CRESLT) 

V    =    DIMAG(CRESLTJ 

RETURN 

END 


AL   «=    1.25D-1 
T  «  8.D+0 
AA   =   0.3 
Xp   =   -4.99 

yF  ■  2.3 

MMMMM   =   24 

LOSC   =    5 

SUM  =   -1.38634D-1 

SUMA   =    3.92274D+0 

TIME   =   1.2126D+0 

ANALYTICAL   =    -1.38634D-1 
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Case  92. 


£ 


-1 


;  [_bi^^.] 


ZI      f  H>0 


[s*  ^«lJ^ 


This  is  known  as  the  Struve  function, 
page  496. 


See  ref.  (1), 


(X,Y,T,AL.3  =  ,iJrV.V"MMMrNCALL  t^CALL) 
RY     *12    °*GE    264    PTPERTS    AND    KAUFMAN 


SU8S3UTIME    VALUE 

THIS     IS    TAPLE    £KI 

IMPLICIT   5 £AL*8    (  a,3 ,0-H,0-Z) 

IMPLIC  IT   C3«PLEX*16    (C  J 

OIMENMCN    NCALU5),     MCALL(5) 

HMMMM    =     MMMMM    *■     1 

ZE5T    =    0.0+0 

CS1    =    0CMPLX(X,Y) 

CS2   =    0CMPLX(AL,ZER3) 

XPTR   =    X*T 

XPTI     =    Y*~ 

CXPT    =    DCMPLX(XPTRtXPTI) 

CS3    =    UCS1-CS1)     *    <CS2*CS2)) 

CALL    CPOWEF.     (CS3,  CQEMCM  ,  5.0-1  , 1  ,  NCALL  ) 

CS4   =     (CDENC1**    ♦    ; 52)  /C  SI 

CALL    CLrG     (C?4tC?5tltMCALLJ 

CNUM    =    CS?*':OEXP(CXPT) 

CRESLT    =    CNUW/COENOM 

U    =    DREAL(C'eSL-) 

V  *    OIMAG(CRESLT) 

RETJUM 

ENO 


AL  =   1.25D-1 
T  ■  8.D+0 
AA  =   0.3 
Xy  »  -4.54 

yF  *  2.2 

MMMMM   =   23 

LOSC   =   k 

SUM  =   8.932^^0-1 

SUMA  =  9.27960D-1 

TIME  =  1.1341D-1 

ANALYTICAL  =  8.932^D-1 
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Case  93 


X 


-1 


Heaviside's  unit  function  f(t  —  a) 


1- 


SUBRDUT 
THIS    IS 

impli:i 

IMPLICI 
MMMMM    = 
CS  =    DC 
XPTR    = 
XPTI     = 
CXPT    = 
CRESLT 
U   *   ORE 
V    =   DIM 
RETJRN 
END 


INE    VALUE     < X,Y,T,AL,U, V, MMMMM) 

TARLE    ENTRY    NUMBER    139    of    LAPLACE    TRANSFORMS    BY    SPIEGEL 
T    R6AL*8     ( A,3t0-H,0-Z) 
t    CnwPLEX*16    (  C  ) 

MMMMM    +    1 
MPLX(X,Y) 
X*(T-AL) 
Y*(T-ALJ 

DCMPLX(XPTR,  XPTI  ) 
=    CDEXPtCXPT  )/CS 
ALCCREiLTJ 
AG(CRESLT) 


AL    =    1.25D-1 
T   =   8.0625D+0 
AA   =   0.2 
Xp    =   -*.6* 

yF  ■  2.2 

MMMMM  =  23 

LOSC  =  5 

SUM  =  1.00000D+0 

SUMA   =   1.U1682D+0 

TIME   =   6.055^0-2 

ANALYTICAL   =    1.00000D+0 
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Case  94 


£ 


-x 


.-* 


=    J.[.*«] 


SUBROUTINE    VALUE     (X ,Y ,T , AL.U, V , VMMMM ) 

THIS    TRANSFORM    WAS    TAKEN    FROM    ROBERTS    AND    KAUFMAN 

IMPLICIT    °EAL*8     (A-H,  O-Z) 

MKMMM    =     MMMMM    ♦     1 

TEM1    =    -l.0+0*4L*AL*X 

TEM2    =     AL#AL*Y 

TEM3    =     4.D+C*(  (  X*X)     +     (Y*Y)J 

TEM4   =     TEM1/TEM3 

TEM5    -    TEM2  /TEM3 

TEM6   =    DEXP(TEM4)*0C0S<TEM5) 

TEM7    =    0EXO(TEM4)*0SIN(TEM5) 

TEM8    =     (TEM6*X)     +     l^^'Y] 

TEM9    =    (TE*I7*X)     -     (TEM6*Y) 

TEM13     =     (X*X)     *     (Y*Y) 

TEM11    =    TEM8/TEM10 

TEM12    =    TEM9/TE^10 

TEM13    =    OEX?(X*TJ*OCnS(Y*T) 

TEM14    =    OEXP{ X*T) *OSIN( Y*T) 

U=     (TEM13*TEMH)    -    <  TEM14*T=M12) 

V    =    (TEMl^*TEMH)    +    (TEMi3*rE"12) 

RETURN 

END 


AL   =   2.50-1 

T   «    1.6D+1 
AA   =   0.2 

xp  =  -3.04 
yF  e  1.8 

MMMMM  =   19 

LOSC    =  8 

SUM  =  7.65198D-1 

SUMA  =  3.98573D+0 

TIME  =  6.36^8D-2 

ANALYTICAL   =   7.65198D-I 
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Case  95 


X 


-i 


i-  fate 


-at 


SUBROUTINE     VALUE    (  X  ,Y  ,T  , AL,U  ,  V  ,  "mmvm) 

THIS    TRANSFORM    WAS    TAKEN    FRO*    ROBERTS    AND    KAUFMAN 

IMPLICIT    REAL*3     I  A, 3, Q-Ht O-ZJ 

IMPLICI  T   COMPLEX*  16    (C) 

MMMMM    =    Mmmmm    +    i 

ZERO    =    O.D+O 

CS1    =    OCMPLXJ  X,  Y) 

CS2    =    DCMPLX  (AL.ZERO) 

CS3   =    CS1    -    CS2 

CS4   =    CS1    +    CS2 

XPTP    =    X*T 

XPTI     =    Y*T 

CXPT    =     DCMPLX(XPTR,XPTI  ) 

CNUM    =    CS3*CS3*CJ  EXP  (CXPT  J 

COENOM    =    CS  l*CSVCS^ 

CRESL"     =    CIUM/CJENOM 

U    =   ORE  AL(CRESLT) 

V  =    OIMAG(CRESLT) 

RETJRN 

ENO 


AL   =    1.25D-1 
T   =   8.D+0 
AA   =    0.2 

yF  ■  2.2 

MMMMM   »   23 

LOSC   =  5 

SUM   =    -4.71518D-1 

SUMA   =   5.0686^D-1 
TIME   =   6.9186D-2 
ANALYTICAL   =    -U.71518D-1 
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